commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From pste...@apache.org
Subject cvs commit: jakarta-commons/math/src/test/org/apache/commons/math/distribution NormalDistributionTest.java
Date Mon, 26 Jan 2004 03:04:31 GMT
psteitz     2004/01/25 19:04:31

  Modified:    math/src/java/org/apache/commons/math/distribution
                        DistributionFactory.java
                        DistributionFactoryImpl.java
  Added:       math/src/java/org/apache/commons/math/distribution
                        NormalCDFAlgorithm.java NormalCDFFastAlgorithm.java
                        NormalCDFPreciseAlgorithm.java
                        NormalDistribution.java NormalDistributionImpl.java
               math/src/test/org/apache/commons/math/distribution
                        NormalDistributionTest.java
  Log:
  Added Normal Distribution implementations and tests contributed by Piotr Kochanski.
  
  Revision  Changes    Path
  1.18      +21 -2     jakarta-commons/math/src/java/org/apache/commons/math/distribution/DistributionFactory.java
  
  Index: DistributionFactory.java
  ===================================================================
  RCS file: /home/cvs/jakarta-commons/math/src/java/org/apache/commons/math/distribution/DistributionFactory.java,v
  retrieving revision 1.17
  retrieving revision 1.18
  diff -u -r1.17 -r1.18
  --- DistributionFactory.java	15 Nov 2003 16:01:35 -0000	1.17
  +++ DistributionFactory.java	26 Jan 2004 03:04:31 -0000	1.18
  @@ -1,7 +1,7 @@
   /* ====================================================================
    * The Apache Software License, Version 1.1
    *
  - * Copyright (c) 2003 The Apache Software Foundation.  All rights
  + * Copyright (c) 2003-2004 The Apache Software Foundation.  All rights
    * reserved.
    *
    * Redistribution and use in source and binary forms, with or without
  @@ -65,6 +65,8 @@
    * <li>Exponential</li>
    * <li>F</li>
    * <li>Gamma</li>
  + * <li>HyperGeometric</li>
  + * <li>Normal</li>
    * <li>Student's t</li>
    * </ul>
    * 
  @@ -164,4 +166,21 @@
       public abstract HypergeometricDistribution
           createHypergeometricDistribution(int populationSize,
               int numberOfSuccesses, int sampleSize);
  + 
  +	/**
  +	 * Create a new normal distribution with the given mean and standard
  +	 * deviation values.
  +	 * @param mean arithmetic mean.
  +	 * @param sd standard deviation.
  +	 * @return a new normal distribution.  
  +	 */           
  +    public abstract NormalDistribution 
  +    	createNormalDistribution(double mean, double sd);
  +    	
  +	/**
  +	 * Create a new normal distribution with the mean equal to zero and standard
  +	 * deviation equal to one.
  +	 * @return a new normal distribution.  
  +	 */               
  +	public abstract NormalDistribution createNormalDistribution();
   }
  
  
  
  1.17      +22 -2     jakarta-commons/math/src/java/org/apache/commons/math/distribution/DistributionFactoryImpl.java
  
  Index: DistributionFactoryImpl.java
  ===================================================================
  RCS file: /home/cvs/jakarta-commons/math/src/java/org/apache/commons/math/distribution/DistributionFactoryImpl.java,v
  retrieving revision 1.16
  retrieving revision 1.17
  diff -u -r1.16 -r1.17
  --- DistributionFactoryImpl.java	19 Nov 2003 03:22:53 -0000	1.16
  +++ DistributionFactoryImpl.java	26 Jan 2004 03:04:31 -0000	1.17
  @@ -1,7 +1,7 @@
   /* ====================================================================
    * The Apache Software License, Version 1.1
    *
  - * Copyright (c) 2003 The Apache Software Foundation.  All rights
  + * Copyright (c) 2003-2004 The Apache Software Foundation.  All rights
    * reserved.
    *
    * Redistribution and use in source and binary forms, with or without
  @@ -153,5 +153,25 @@
           return new HypergeometricDistributionImpl(populationSize,
               numberOfSuccesses, sampleSize);
       }
  +
  +	/**
  +	 * Create a new normal distribution with the given mean and standard
  +	 * deviation values.
  +	 * @param mean arithmetic mean.
  +	 * @param sd standard deviation.
  +	 * @return a new normal distribution.  
  +	 */   
  +	public NormalDistribution createNormalDistribution(double mean, double sd) {
  +		return new NormalDistributionImpl(mean, sd);
  +	}
  +
  +	/**
  +	 * Create a new normal distribution with the mean equal to zero and standard
  +	 * deviation equal to one.
  +	 * @return a new normal distribution.  
  +	 */ 
  +	public NormalDistribution createNormalDistribution() {
  +		return new NormalDistributionImpl();
  +	}
   
   }
  
  
  
  1.1                  jakarta-commons/math/src/java/org/apache/commons/math/distribution/NormalCDFAlgorithm.java
  
  Index: NormalCDFAlgorithm.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2004 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowledgement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowledgement may appear in the software itself,
   *    if and wherever such third-party acknowledgements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their name without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  
  package org.apache.commons.math.distribution;
  
  /**
   * Interface for various algorithms, that calculate CDF 
   * (cumulative distribution function).
   */
  interface NormalCDFAlgorithm {
  	/**
  	 * For normal distribution X, this method returns P(X &lt; x).
  	 * @param x the value at which the CDF is evaluated.
  	 * @return CDF for this distribution. 
  	 */	
  	double cdf(double x);
  }
  
  
  
  1.1                  jakarta-commons/math/src/java/org/apache/commons/math/distribution/NormalCDFFastAlgorithm.java
  
  Index: NormalCDFFastAlgorithm.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2004 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowledgement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowledgement may appear in the software itself,
   *    if and wherever such third-party acknowledgements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their name without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  
  package org.apache.commons.math.distribution;
  
  /**
   * Implements Hill algorithm to calculate CDF (cumulative distribution function)
   * for Normal (Gauss) distribution.<p>
   * References:<p> 
   * I. D. Hill, Applied Statistics, 
   * <a href="http://lib.stat.cmu.edu/apstat/66">Algorithm AS 66: 
   * The Normal Integral</a>. 
   * <p>This algorithm is faster <b>but</b> less precise then the one implemented
   * in {@link org.apache.commons.math.distribution.NormalCDFPreciseAlgorithm}. 
   */
  class NormalCDFFastAlgorithm implements NormalCDFAlgorithm {
  
  	/**
  	 * Calculates P(X &lt; <code>x</code>).
  	 * @param x the value at which the CDF is evaluated.
  	 * @return CDF evaluted at <code>x</code>. 
  	 */
  	public double cdf(double x) {
  		final double ltone = 7.0,
  				utzero = 18.66,
  				con = 1.28,
  				p = 0.398942280444,
  				q = 0.399903438504,
  				r = 0.398942280385,
  				a1 = 5.75885480458,
  				a2 = 2.62433121679,
  				a3 = 5.92885724438,
  				b1 = -29.8213557808,
  				b2 = 48.6959930692,
  				c1 = -3.8052e-8,
  				c2 = 3.98064794e-4,
  				c3 = -0.151679116635,
  				c4 = 4.8385912808,
  				c5 = 0.742380924027,
  				c6 = 3.99019417011,
  				d1 = 1.00000615302,
  				d2 = 1.98615381364,
  				d3 = 5.29330324926,
  				d4 = -15.1508972451,
  				d5 = 30.789933034;
  		boolean upper = false; //kept to preserve similarity with the original 
  							   //algorithm
  		double y, alnorm;
  
  		if (x < 0) {
  			upper = !upper;
  			x = -x;
  		}
  		if (x <= ltone || upper && x <= utzero) {
  			y = 0.5 * x * x;
  			if(x>con) {
  				alnorm=r*Math.exp(-y)/(x+c1+d1/(x+c2+d2/(x+c3+d3/(x+
  						c4+d4/(x+c5+d5/(x+c6))))));
  			}else {
  				alnorm=0.5-x*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3))));
  			}
  		} else {
  			alnorm = 0;
  		}
  		if (!upper)
  			alnorm = 1 - alnorm;
  		return alnorm;
  	}
  }
  
  
  
  1.1                  jakarta-commons/math/src/java/org/apache/commons/math/distribution/NormalCDFPreciseAlgorithm.java
  
  Index: NormalCDFPreciseAlgorithm.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2004 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowledgement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowledgement may appear in the software itself,
   *    if and wherever such third-party acknowledgements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their name without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  
  package org.apache.commons.math.distribution;
  
  /**
   * Implements Cody algorithm to calculate CDF (cumulative distribution function)
   * for Normal (Gauss) distribution.<p> 
   * Provided implementation is adapted from 
   * <a href="http://www.r-project.org/">R statistical package</a> function
   * <code>pnorm(...)</code>.<p>
   * References:
   * <ul> 
   * <li>Cody's algorithm: Cody, W. D. (1993).<br>
   *	<a href="http://www.acm.org/toms/cgi-bin/TOMSbibget?Cody:1993:ASE">
   *  ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of
   *	Special Function Routines and Test Drivers".
   *	ACM Transactions on Mathematical Software. 19, 22-32</a>.</li>
   * <li>FORTRAN code is available at 
   * <a href="http://www.netlib.org/toms/715">http://www.netlib.org/toms/715</a>
   * </li>
   * </ul>
   */
  class NormalCDFPreciseAlgorithm implements NormalCDFAlgorithm{
  	/**
  	 * Calculates P(X &lt; <code>x</code>).
  	 * @param x the value at which the CDF is evaluated.
  	 * @return CDF evaluted at <code>x</code>. 
  	 */	
  	public double cdf(double x){
  		final double[] a = {
  			2.2352520354606839287, 1.6102823106855587881e2,
  			1.0676894854603709582e3, 1.8154981253343561249e4,
  			6.5682337918207449113e-2
  		};
  		final double[] b = {
  			4.7202581904688241870e1, 9.7609855173777669322e2,
  			1.0260932208618978205e4, 4.5507789335026729956e4
  		};
  		final double[] c = {
  			3.9894151208813466764e-1, 8.8831497943883759412,
  			9.3506656132177855979e1, 5.9727027639480026226e2,
  			2.4945375852903726711e3,6.8481904505362823326e3,
  			1.1602651437647350124e4, 9.8427148383839780218e3,
  			1.0765576773720192317e-8
  		};
  		final double[] d = {
  			2.2266688044328115691e1, 2.3538790178262499861e2,
  			1.5193775994075548050e3, 6.4855582982667607550e3, 
  			1.8615571640885098091e4, 3.4900952721145977266e4,
  			3.8912003286093271411e4, 1.9685429676859990727e4
  		};
  		final double[] p = {
  			2.1589853405795699e-1, 1.274011611602473639e-1,
  			2.2235277870649807e-2, 1.421619193227893466e-3,
  			2.9112874951168792e-5, 2.307344176494017303e-2
  		};
  		final double[] q = {
  			1.28426009614491121, 4.68238212480865118e-1,
  			6.59881378689285515e-2, 3.78239633202758244e-3,
  			7.29751555083966205e-5
  		};
  		
  		final double sixten = 16.0;
  		// 1/sqrt(2pi)
  		final double M_1_SQRT_2PI = 0.398942280401432677939946059934;
  		// sqrt(32)	
  		final double M_SQRT_32 = 5.656854249492380195206754896838;
  
  		final double eps = 0.5*Math.pow(2,(double)(1-60));
  		final double min = Double.MIN_VALUE;
  		
  		double ccum = 0, cum = 0;
  		double del, temp, xden = 0, xnum = 0, xsq = 0;
  		
  		int i;
  
  		double y = Math.abs(x);
  		//1st case: |x| <= qnorm(3/4)
  		if(y <= 0.67448975) { 
  			if(y > eps) xsq = x*x;
  			xnum = a[4]*xsq;
  			xden = xsq;
  			for(i=0; i < 3; i++) {
  				xnum = (xnum + a[i])*xsq;
  				xden = (xden + b[i])*xsq;
  			}
  			temp = x*(xnum + a[3])/(xden + b[3]);
  			cum = 0.5 + temp;
  			ccum = 0.5 - temp;
  		}
  		//2nd case: qnorm(3/4)1 <= |x| <= sqrt(32)
  		else if(y <= M_SQRT_32) {
  			xnum = c[8]*y;
  			xden = y;
  			for(i=0; i < 7; i++) {
  				xnum = (xnum + c[i])*y;
  				xden = (xden + d[i])*y;
  			}
  			temp = (xnum + c[7])/(xden + d[7]);
  			xsq = Math.floor(y*sixten)/sixten;
  			del = (y - xsq)*(y + xsq);
  			cum = Math.exp(-(xsq*xsq*0.5))*Math.exp(-(del*0.5))*temp;
  			ccum = 1.0 - cum;
  			if(x > 0.0) {
  				temp = cum;
  				cum = ccum;
  				ccum = temp;
  			}
  		}
  		//3rd case: -37.5193 < x  &&  x < 8.2924 || -8.2924  < x  &&
 x < 37.5193
  		else if(-37.5193 < x  &&  x < 8.2924 || -8.2924  < x  &&  x <
37.5193) {
  			xsq = 1.0/(x*x);
  			xnum = p[5]*xsq;
  			xden = xsq;
  			for(i=0; i < 4; i++) {
  				xnum = (xnum+p[i])*xsq;
  				xden = (xden+q[i])*xsq;
  			}
  			temp = xsq*(xnum + p[4])/(xden + q[4]);
  			temp = (M_1_SQRT_2PI - temp)/y;
  			xsq = Math.floor(x*sixten)/sixten;
  			del = (x - xsq)*(x + xsq);
  			cum = Math.exp(-(xsq*xsq*0.5))*Math.exp(-(del*0.5))*temp;
  			ccum = 1.0 - cum;
  			if(x > 0.0) {
  				temp = cum;
  				cum = ccum;
  				ccum = cum;
  			}
  		//4th case: high x values	
  		}else{
  			if(x > 0) {	
  				cum = 1.0; 
  				ccum = 0.0;	
  			}else {
  				cum = 0.0; 
  				ccum = 1.0;	
  			}
  		}
  		if(cum < min) cum = 0.0;
  		if(ccum < min) ccum = 0.0;
  		
  		return cum;
  	}
  }
  
  
  
  1.1                  jakarta-commons/math/src/java/org/apache/commons/math/distribution/NormalDistribution.java
  
  Index: NormalDistribution.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2004 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowledgement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowledgement may appear in the software itself,
   *    if and wherever such third-party acknowledgements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their name without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  
  package org.apache.commons.math.distribution;
  
  /**
   * Normal (Gauss) Distribution.
   * Instances of NormalDistribution objects should be created using
   * {@link DistributionFactory#createNormalDistribution(double, double)}.<p>
   * 
   * References:<p>
   * <ul>
   * <li><a href="http://mathworld.wolfram.com/NormalDistribution.html">
   * Normal Distribution</a></li>
   * </ul>
   * 
   */
  public interface NormalDistribution extends ContinuousDistribution {
  	/**
  	 * Access the mean.
  	 * @return mean for this distribution
  	 */
  	double getMean();
  	/**
  	 * Modify the mean.
  	 * @param mean for this distribution
  	 */
  	void setMean(double mean);
  	/**
  	 * Access the standard deviation.
  	 * @return standard deviation for this distribution
  	 */
  	double getStandardDeviation();
  	/**
  	 * Modify the standard deviation.
  	 * @param sd standard deviation for this distribution
  	 */
  	void setStandardDeviation(double sd);
  	
  	/**
  	 * Access algorithm used to calculate cummulative probability
  	 * @return cdfAlgorithm the value of cummulative probability
  	 */
  	public NormalCDFAlgorithm getCdfAlgorithm();
  
  	/**
  	 * Modify the algorithm used to calculate cummulative probability
  	 * @param normalCDF the algorithm used to calculate cummulative probability
  	 */
  	public void setCdfAlgorithm(NormalCDFAlgorithm normalCDF);
  }
  
  
  
  1.1                  jakarta-commons/math/src/java/org/apache/commons/math/distribution/NormalDistributionImpl.java
  
  Index: NormalDistributionImpl.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2004 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowledgement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowledgement may appear in the software itself,
   *    if and wherever such third-party acknowledgements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their name without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  
  package org.apache.commons.math.distribution;
  
  import java.io.Serializable;
  
  /**
   * Default implementation of
   * {@link org.apache.commons.math.distribution.NormalDistribution}.<p>
   * You can choose the algorithm used to calculate cummulative probability
   * using method {@link #setCdfAlgorithm}. The deafault is the Cody algorithm 
   * {@link org.apache.commons.math.distribution.NormalCDFPreciseAlgorithm}
   */
  public class NormalDistributionImpl extends AbstractContinuousDistribution 
  		implements NormalDistribution, Serializable {
  	private double mean = 0;
  	private double standardDeviation = 1;
  	private NormalCDFAlgorithm cdfAlgorithm = new NormalCDFPreciseAlgorithm();
  	
  	/**
  	 * Create a normal distribution using the given mean and standard deviation.
  	 * @param mean mean for this distribution
  	 * @param sd standard deviation for this distribution
  	 */
  	public NormalDistributionImpl(double mean, double sd){
  		super();
  		setMean(mean);
  		setStandardDeviation(sd);
  	}
  	/**
  	 * Creates normal distribution with the mean equal to zero and standard
  	 * deviation equal to one. 
  	 */
  	public NormalDistributionImpl(){
  		super();
  		setMean(0.0);
  		setStandardDeviation(1.0);
  	}	
  	/**
  	 * Access the mean.
  	 * @return mean for this distribution
  	 */	
  	public double getMean() {
  		return mean;
  	}
  	/**
  	 * Modify the mean.
  	 * @param mean for this distribution
  	 */
  	public void setMean(double mean) {
  		this.mean = mean;
  	}
  
  	/**
  	 * Access the standard deviation.
  	 * @return standard deviation for this distribution
  	 */
  	public double getStandardDeviation() {
  		return standardDeviation;
  	}
  
  	/**
  	 * Modify the standard deviation.
  	 * @param sd standard deviation for this distribution
  	 */
  	public void setStandardDeviation(double sd) {
  		if (sd < 0.0) {
  			throw new IllegalArgumentException("Standard deviation must be" +
				"positive or zero.");
  		}		
  		standardDeviation = sd;
  	}
  
  	/**
  	 * For this disbution, X, this method returns P(X &lt; <code>x</code>).
  	 * @param x the value at which the CDF is evaluated.
  	 * @return CDF evaluted at <code>x</code>. 
  	 */
  	public double cummulativeProbability(double x) {
  		double z = x;
  		if(standardDeviation > 0){
  			z = (x - mean)/standardDeviation;
  		}else{
  			return 0.0;
  		}
  		return cdfAlgorithm.cdf(z);
  	}
  
  
  	/**
  	 * For this distribution, X, this method returns the critical point x, such
  	 * that P(X &lt; x) = <code>p</code>.<p>
  	 * Provided implementation is adopted from 
       * <a href="http://www.r-project.org/">R statistical package</a> function
       * <code>qnorm(...)</code>.<p>
  	 * References:
  	 * <ul>
  	 * <li>
  	 *  Beasley, J. D. and S. G. Springer (1977).
  	 *  <a href="http://lib.stat.cmu.edu/apstat/111">
  	 *	Algorithm AS 111: The percentage points of the normal distribution</a>,
  	 *	Applied Statistics, 26, 118-121.
  	 * </li>
  	 * <li>
  	 *  Wichura, M.J. (1988).
  	 *  <a href="http://lib.stat.cmu.edu/apstat/241">
  	 *  Algorithm AS 241: The Percentage Points of the Normal Distribution.</a>
  	 *  Applied Statistics, 37, 477-484.
  	 * </li>
  	 * </ul>
  	 *
  	 * @param p the desired probability
  	 * @return x, such that P(X &lt; x) = <code>p</code>
  	 */
  	public double inverseCummulativeProbability(double p) {
  		if (p < 0.0 || p > 1.0) {
  			throw new IllegalArgumentException("p must be between 0.0 and 1.0, inclusive.");
  		}
  		
  		//TODO is this ok?
  		if(standardDeviation == 0){
  			return mean;
  		}
  		
  		double r, val;		
  		double q = p - 0.5;
  
  		if (Math.abs(q) <= .425) {/* 0.075 <= p <= 0.925 */
  			r = 0.180625 - q*q;
  			val =
  				q * (((((((r * 2509.0809287301226727 +
  						   33430.575583588128105) * r + 67265.770927008700853) * r +
  						 45921.953931549871457) * r + 13731.693765509461125) * r +
  					   1971.5909503065514427) * r + 133.14166789178437745) * r +
  					 3.387132872796366608)
  				/ (((((((r * 5226.495278852854561 +
  						 28729.085735721942674) * r + 39307.89580009271061) * r +
  					   21213.794301586595867) * r + 5394.1960214247511077) * r +
  					 687.1870074920579083) * r + 42.313330701600911252) * r + 1.);
  		}else { //closer than 0.075 from {0,1} boundary
  		if (q > 0)
  			r = 1 - p;
  		else
  			r = p;
  		r = Math.sqrt(- Math.log(r));
  		if (r <= 5.0) {
  			r += -1.6;
  			val = (((((((r * 7.7454501427834140764e-4 +
  					   0.0227238449892691845833) * r + 0.24178072517745061177) *
  					 r + 1.27045825245236838258) * r +
  					3.64784832476320460504) * r + 5.7694972214606914055) *
  				  r + 4.6303378461565452959) * r +
  				 1.42343711074968357734)
  				/ (((((((r *
  						 1.05075007164441684324e-9 + 5.475938084995344946e-4) *
  						r + 0.0151986665636164571966) * r +
  					   0.14810397642748007459) * r + 0.68976733498510000455) *
  					 r + 1.6763848301838038494) * r +
  					2.05319162663775882187) * r + 1.0);
  		}else { //very close to  0 or 1
  			r += -5.;
  			val = (((((((r * 2.01033439929228813265e-7 +
  					   2.71155556874348757815e-5) * r +
  					  0.0012426609473880784386) * r + 0.026532189526576123093) *
  					r + 0.29656057182850489123) * r +
  				   1.7848265399172913358) * r + 5.4637849111641143699) *
  				 r + 6.6579046435011037772)
  				/ (((((((r *
  						 2.04426310338993978564e-15 + 1.4215117583164458887e-7)*
  						r + 1.8463183175100546818e-5) * r +
  					   7.868691311456132591e-4) * r + 0.0148753612908506148525)
  					 * r + 0.13692988092273580531) * r +
  					0.59983220655588793769) * r + 1.0);
  		}
  		if(q < 0.0)
  			val = -val;
  		}
  		return mean + standardDeviation*val;
  	}
  
  
  	/**
  	 * Access algorithm used to calculate cummulative probability
  	 * @return cdfAlgorithm the value of cummulative probability
  	 */
  	public NormalCDFAlgorithm getCdfAlgorithm() {
  		return cdfAlgorithm;
  	}
  
  
  	/**
  	 * Modify the algorithm used to calculate cummulative probability
  	 * @param normalCDF the algorithm used to calculate cummulative probability
  	 */
  	public void setCdfAlgorithm(NormalCDFAlgorithm normalCDF) {
  		cdfAlgorithm = normalCDF;
  	}
  
  	
  	/**
  	 * Access the domain value lower bound, based on <code>p</code>, used to
  	 * bracket a CDF root.  This method is used by
  	 * {@link #inverseCummulativeProbability(double)} to find critical values.
  	 * 
  	 * @param p the desired probability for the critical value
  	 * @return domain value lower bound, i.e.
  	 *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>

  	 */
  	protected double getDomainLowerBound(double p) {
  		return -Double.MAX_VALUE;
  	}
  
  	/**
  	 * Access the domain value upper bound, based on <code>p</code>, used to
  	 * bracket a CDF root.  This method is used by
  	 * {@link #inverseCummulativeProbability(double)} to find critical values.
  	 * 
  	 * @param p the desired probability for the critical value
  	 * @return domain value upper bound, i.e.
  	 *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>

  	 */
  	protected double getDomainUpperBound(double p) {
  		return Double.MAX_VALUE;
  	}
  
  	/**
  	 * Access the initial domain value, based on <code>p</code>, used to
  	 * bracket a CDF root.  This method is used by
  	 * {@link #inverseCummulativeProbability(double)} to find critical values.
  	 * 
  	 * @param p the desired probability for the critical value
  	 * @return initial domain value
  	 */
  	protected double getInitialDomain(double p) {
  		return 0.0;
  	}
  
  
  }
  
  
  
  1.1                  jakarta-commons/math/src/test/org/apache/commons/math/distribution/NormalDistributionTest.java
  
  Index: NormalDistributionTest.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2004 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   * 
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowledgement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowledgement may appear in the software itself,
   *    if and wherever such third-party acknowledgements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their name without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  
  package org.apache.commons.math.distribution;
  
  import org.apache.commons.math.MathException;
  
  import junit.framework.TestCase;
  
  /**
   *  Tests for NormalDistribution implementation
   * 
   * "True" results are taken from R - the same as in Mathematica
   *
   */
  public class NormalDistributionTest extends TestCase {
  	
  	private NormalDistribution z;
  	private static final double PRECISION = 10e-6;	
  	private static final double M = 2.1;
  	private static final double SD = 1.4;
  	
  	/**
  	 * Constructor for NormalDistributionTest.
  	 * @param arg0
  	 */
  	public NormalDistributionTest(String arg0) {
  		super(arg0);
  	}
  
  	public static void main(String[] args) {
  		junit.swingui.TestRunner.run(NormalDistributionTest.class);
  	}
  
  	protected void setUp() throws Exception {
  		super.setUp();
  		z = DistributionFactory.newInstance().createNormalDistribution(M, SD);
  	}
  
  	protected void tearDown() throws Exception {
  		super.tearDown();
  		z = null;
  	}
  
  	public void testCummulativeProbabilitydoubleM_MINUS_2SD() throws MathException {
  		testProbability(M - 2*SD, 0.02275013);
  	}
  
  	public void testCummulativeProbabilitydoubleM_MINUS_SD() throws MathException {
  		testProbability(M - SD, 0.1586553);
  	}
  
  	public void testCummulativeProbabilitydoubleM() throws MathException {
  		testProbability(M, 0.5);
  	}
  
  	public void testCummulativeProbabilitydoubleM_PLUS_SD() throws MathException {
  		testProbability(M + SD, 0.8413447);
  	}
  	
  	public void testCummulativeProbabilitydoubleM_PLUS_2SD() throws MathException {
  		testProbability(M + 2*SD, 0.9772499);
  	}
  	
  	public void testCummulativeProbabilitydoubleM_PLUS_3SD() throws MathException {
  		testProbability(M + 3*SD, 0.9986501);
  	}
  	
  	public void testCummulativeProbabilitydoubleM_PLUS_4SD() throws MathException {
  		testProbability(M + 4*SD, 0.9999683);
  	}
  	
  	public void testCummulativeProbabilitydoubleM_PLUS_5SD() throws MathException {
  		testProbability(M + 5*SD, 0.9999997);
  	}
  	
  	public void testInverseCummulativeProbability0() throws MathException {
  		assertEquals(Double.isNaN(z.inverseCummulativeProbability(0.0)), true);
  	}
  
  	public void testInverseCummulativeProbability001() throws MathException {
  		testValue(-2.226325, .001);
  	}
  
  	public void testInverseCumulativeProbability010() throws MathException{
  		testValue(-1.156887, .010);
  	}
  
  	public void testInverseCumulativeProbability025() throws MathException{
  		testValue(-0.6439496, .025);
  	}
  
  	public void testInverseCumulativeProbability050() throws MathException{
  		testValue(-0.2027951, .050);
  	}
  
  	public void testInverseCumulativeProbability100() throws MathException{
  		testValue(0.3058278, .100);
  	}
  
  	public void testInverseCumulativeProbability900() throws MathException{
  		testValue(3.894172, .900);
  	}
  
  	public void testInverseCumulativeProbability950() throws MathException{
  		testValue(4.402795, .950);
  	}
  
  	public void testInverseCumulativeProbability975() throws MathException{
  		testValue(4.84395, .975);
  	}
  
  	public void testInverseCumulativeProbability990() throws MathException{
  		testValue(5.356887, .990);
  	}
  
  	public void testInverseCummulativeProbability999() throws MathException{
  		testValue(6.426325, .999);
  	}
  
  	public void testInverseCummulativeProbability1() throws MathException {
  		assertEquals(Double.isNaN(z.inverseCummulativeProbability(1.0)), true);
  	}
  
  	public void testGetMean() {
  		assertEquals(M, z.getMean(), 0);
  	}
  
  	public void testSetMean() throws MathException {
  		double mu = Math.random();
  		z.setMean(mu);
  		assertEquals(mu, z.getMean(), 0);
  		assertEquals(0.5d, z.cummulativeProbability(mu), PRECISION);
  	}
  
  	public void testGetStandardDeviation() {
  		assertEquals(SD, z.getStandardDeviation(), 0);	
  	}
  
  	public void testSetStandardDeviation() throws MathException{
  		double sigma = 0.1d + Math.random();
  		z.setStandardDeviation(sigma);
  		assertEquals(sigma, z.getStandardDeviation(), 0);
  		assertEquals(0.84134475, z.cummulativeProbability(z.getMean() + z.getStandardDeviation()),
PRECISION );
  	}
  
  	public void testGetCdfAlgorithm() {
  		assertTrue(z.getCdfAlgorithm() != null);
  	}
  
  	public void testSetCdfAlgorithm() {
  		z.setCdfAlgorithm(new NormalCDFFastAlgorithm());
  		assertTrue(z.getCdfAlgorithm() instanceof NormalCDFFastAlgorithm);
  	}
  	
  	private void testProbability(double x, double expected) throws MathException {
  		double actual = Double.NaN;
  		z.setCdfAlgorithm(new NormalCDFPreciseAlgorithm());
  		actual =  z.cummulativeProbability(x);
  		assertEquals(expected, actual, PRECISION);
  		z.setCdfAlgorithm(new NormalCDFFastAlgorithm());
  		actual =  z.cummulativeProbability(x);
  		assertEquals(expected, actual, PRECISION);
  	}
  
  	private void testValue(double expected, double p) throws MathException {
  		double actual = z.inverseCummulativeProbability(p);
  		assertEquals(expected, actual, PRECISION);
  	}
  
  }
  
  
  

---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message