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 [EMAIL PROTECTED]
*
* 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 < 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 [EMAIL PROTECTED]
*
* 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 [EMAIL PROTECTED]
org.apache.commons.math.distribution.NormalCDFPreciseAlgorithm}.
*/
class NormalCDFFastAlgorithm implements NormalCDFAlgorithm {
/**
* Calculates P(X < <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 [EMAIL PROTECTED]
*
* 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 < <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 [EMAIL PROTECTED]
*
* 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
* [EMAIL PROTECTED] 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 [EMAIL PROTECTED]
*
* 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
* [EMAIL PROTECTED] org.apache.commons.math.distribution.NormalDistribution}.<p>
* You can choose the algorithm used to calculate cummulative probability
* using method [EMAIL PROTECTED] #setCdfAlgorithm}. The deafault is the Cody
algorithm
* [EMAIL PROTECTED] 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 < <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 < 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 < 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
* [EMAIL PROTECTED] #inverseCummulativeProbability(double)} to find critical
values.
*
* @param p the desired probability for the critical value
* @return domain value lower bound, i.e.
* P(X < <i>lower bound</i>) < <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
* [EMAIL PROTECTED] #inverseCummulativeProbability(double)} to find critical
values.
*
* @param p the desired probability for the critical value
* @return domain value upper bound, i.e.
* P(X < <i>upper bound</i>) > <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
* [EMAIL PROTECTED] #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 [EMAIL PROTECTED]
*
* 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: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
