Added Normal Distribution implementations and tests contributed by Piotr Kochanski.

git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@141065 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Phil Steitz 2004-01-26 03:04:31 +00:00
parent 1b96f28e41
commit fcde02ecce
5 changed files with 657 additions and 4 deletions

View File

@ -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 @@ import org.apache.commons.discovery.tools.DiscoverClass;
* <li>Exponential</li>
* <li>F</li>
* <li>Gamma</li>
* <li>HyperGeometric</li>
* <li>Normal</li>
* <li>Student's t</li>
* </ul>
*
@ -75,7 +77,7 @@ import org.apache.commons.discovery.tools.DiscoverClass;
* ChiSquaredDistribution chi = factory.createChiSquareDistribution(5.0);
* </pre>
*
* @version $Revision: 1.17 $ $Date: 2003/11/15 16:01:35 $
* @version $Revision: 1.18 $ $Date: 2004/01/26 03:04:31 $
*/
public abstract class DistributionFactory {
/**
@ -164,4 +166,21 @@ public abstract class DistributionFactory {
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();
}

View File

@ -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
@ -60,7 +60,7 @@ import java.io.Serializable;
* A concrete distribution factory. This is the default factory used by
* Commons-Math.
*
* @version $Revision: 1.16 $ $Date: 2003/11/19 03:22:53 $
* @version $Revision: 1.17 $ $Date: 2004/01/26 03:04:31 $
*/
public class DistributionFactoryImpl extends DistributionFactory implements Serializable {
/**
@ -154,4 +154,24 @@ public class DistributionFactoryImpl extends DistributionFactory implements Seri
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();
}
}

View File

@ -0,0 +1,102 @@
/* ====================================================================
* 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);
}

View File

@ -0,0 +1,290 @@
/* ====================================================================
* 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;
}
}

View File

@ -0,0 +1,222 @@
/* ====================================================================
* 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);
}
}