Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

PR: http://nagoya.apache.org/bugzilla/show_bug.cgi?id=20766 

Submitted by:	brent@worden.org


git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@140905 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Mark R. Diggory 2003-06-14 04:17:49 +00:00
parent 3a7fab66ae
commit d20ee8ab0e
5 changed files with 364 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
src/java/org/apache/commons/math/distribution/TDistributionImpl.java
View File

@ -103,18 +103,21 @@ public class TDistributionImpl
* @return CDF evaluted at <code>x</code>.
*/
public double cummulativeProbability(double x) {
double t = Beta.regularizedBeta(
getDegreesOfFreedom() / (getDegreesOfFreedom() + (x * x)),
0.5 * getDegreesOfFreedom(), 0.5);
double ret;
if(x < 0.0){
ret = 0.5 * t;
} else if(x > 0.0){
ret = 1.0 - 0.5 * t;
} else {
if(x == 0.0){
ret = 0.5;
} else {
double t = Beta.regularizedBeta(
getDegreesOfFreedom() / (getDegreesOfFreedom() + (x * x)),
0.5 * getDegreesOfFreedom(), 0.5);
if(x < 0.0){
ret = 0.5 * t;
} else {
ret = 1.0 - 0.5 * t;
}
}
return ret;
}

20
src/java/org/apache/commons/math/special/Beta.java
View File

@ -139,17 +139,9 @@ public class Beta {
public static double regularizedBeta(double x, final double a, final double b, double epsilon, int maxIterations) {
double ret;
if (a <= 0.0) {
throw new IllegalArgumentException("a must be positive");
} else if (b <= 0.0) {
throw new IllegalArgumentException("b must be positive");
} else if (x < 0.0 || x > 1.0) {
throw new IllegalArgumentException(
"x must be between 0.0 and 1.0, inclusive");
} else if(x == 0.0){
ret = 0.0;
} else if(x == 1.0){
ret = 1.0;
if (Double.isNaN(x) || Double.isNaN(a) || Double.isNaN(b) || (x < 0)
|| (x > 1) || (a <= 0.0) || (b <= 0.0)) {
ret = Double.NaN;
} else {
ContinuedFraction fraction = new ContinuedFraction() {
protected double getB(int n, double x) {
@ -228,10 +220,8 @@ public class Beta {
public static double logBeta(double a, double b, double epsilon, int maxIterations) {
double ret;
if (a <= 0.0) {
throw new IllegalArgumentException("a must be positive");
} else if (b <= 0.0) {
throw new IllegalArgumentException("b must be positive");
if (Double.isNaN(a) || Double.isNaN(b) || (a <= 0.0) || (b <= 0.0)) {
ret = Double.NaN;
} else {
ret = Gamma.logGamma(a, epsilon, maxIterations) + Gamma.logGamma(b, epsilon, maxIterations)
- Gamma.logGamma(a + b, epsilon, maxIterations);

214
src/java/org/apache/commons/math/stat/UnivariateImpl.java
View File

@ -70,7 +70,7 @@ import org.apache.commons.math.FixedDoubleArray;
* @author <a href="mailto:tobrien@apache.org">Tim O'Brien</a>
* @author <a href="mailto:mdiggory@apache.org">Mark Diggory</a>
* @author Brent Worden
* @version $Revision: 1.3 $ $Date: 2003/06/04 12:23:44 $
* @version $Revision: 1.4 $ $Date: 2003/06/14 04:17:49 $
*
*/
public class UnivariateImpl implements Univariate, Serializable {
@ -89,11 +89,11 @@ public class UnivariateImpl implements Univariate, Serializable {
/** running sum of squares that have been added */
private double sumsq = 0.0;
/** running sum of 3rd powers that have been added */
private double sumCube = 0.0;
/** running sum of 3rd powers that have been added */
private double sumCube = 0.0;
/** running sum of 4th powers that have been added */
private double sumQuad = 0.0;
/** running sum of 4th powers that have been added */
private double sumQuad = 0.0;
/** count of values that have been added */
private int n = 0;
@ -120,17 +120,17 @@ public class UnivariateImpl implements Univariate, Serializable {
/**
* @see org.apache.commons.math.stat.Univariate#addValue(double)
*/
public void addValue(double v) {
* @see org.apache.commons.math.stat.Univariate#addValue(double)
*/
public void addValue(double v) {
insertValue(v);
}
/**
* @see org.apache.commons.math.stat.Univariate#getMean()
*/
public double getMean() {
* @see org.apache.commons.math.stat.Univariate#getMean()
*/
public double getMean() {
if (n == 0) {
return Double.NaN;
} else {
@ -140,97 +140,97 @@ public class UnivariateImpl implements Univariate, Serializable {
/**
* @see org.apache.commons.math.stat.Univariate#getGeometricMean()
*/
public double getGeometricMean() {
* @see org.apache.commons.math.stat.Univariate#getGeometricMean()
*/
public double getGeometricMean() {
if ((product <= 0.0) || (n == 0)) {
return Double.NaN;
} else {
return Math.pow(product,( 1.0/(double)n ) );
return Math.pow(product,( 1.0 / (double) n ) );
}
}
/**
* @see org.apache.commons.math.stat.Univariate#getProduct()
*/
public double getProduct() {
* @see org.apache.commons.math.stat.Univariate#getProduct()
*/
public double getProduct() {
return product;
}
/**
* @see org.apache.commons.math.stat.Univariate#getStandardDeviation()
*/
public double getStandardDeviation() {
double variance = getVariance();
if ((variance == 0.0) || (variance == Double.NaN)) {
return variance;
} else {
return Math.sqrt(variance);
}
}
/**
* @see org.apache.commons.math.stat.Univariate#getStandardDeviation()
*/
public double getStandardDeviation() {
double variance = getVariance();
if ((variance == 0.0) || (variance == Double.NaN)) {
return variance;
} else {
return Math.sqrt(variance);
}
}
/**
* Returns the variance of the values that have been added as described by
* <a href="http://mathworld.wolfram.com/k-Statistic.html">Equation (5) for k-Statistics</a>.
*
* @return The variance of a set of values. Double.NaN is returned for
* an empty set of values and 0.0 is returned for a &lt;= 1 value set.
*/
public double getVariance() {
double variance = Double.NaN;
/**
* Returns the variance of the values that have been added as described by
* <a href="http://mathworld.wolfram.com/k-Statistic.html">Equation (5) for k-Statistics</a>.
*
* @return The variance of a set of values. Double.NaN is returned for
* an empty set of values and 0.0 is returned for a &lt;= 1 value set.
*/
public double getVariance() {
double variance = Double.NaN;
if( n == 1 ) {
variance = 0.0;
} else if( n > 1 ) {
variance = (((double)n)*sumsq - (sum * sum)) / (double) (n * (n - 1));
}
if( n == 1 ) {
variance = 0.0;
} else if( n > 1 ) {
variance = (((double) n) * sumsq - (sum * sum)) / (double) (n * (n - 1));
}
return variance < 0 ? 0.0 : variance;
}
return variance < 0 ? 0.0 : variance;
}
/**
* Returns the skewness of the values that have been added as described by
/**
* Returns the skewness of the values that have been added as described by
* <a href="http://mathworld.wolfram.com/k-Statistic.html">Equation (6) for k-Statistics</a>.
*
* @return The skew of a set of values. Double.NaN is returned for
* an empty set of values and 0.0 is returned for a &lt;= 2 value set.
*/
public double getSkewness() {
* @return The skew of a set of values. Double.NaN is returned for
* an empty set of values and 0.0 is returned for a &lt;= 2 value set.
*/
public double getSkewness() {
if( n < 1) return Double.NaN;
if( n <= 2 ) return 0.0;
if( n < 1) return Double.NaN;
if( n <= 2 ) return 0.0;
return ( 2*Math.pow(sum,3) - 3*sum*sumsq + ((double)(n*n))*sumCube ) /
( (double)(n*(n-1)*(n-2)) ) ;
}
return ( 2 * Math.pow(sum, 3) - 3 * sum * sumsq + ((double) (n * n)) * sumCube ) /
( (double) (n * (n - 1) * (n - 2)) ) ;
}
/**
* Returns the kurtosis of the values that have been added as described by
/**
* Returns the kurtosis of the values that have been added as described by
* <a href="http://mathworld.wolfram.com/k-Statistic.html">Equation (7) for k-Statistics</a>.
*
* @return The kurtosis of a set of values. Double.NaN is returned for
* an empty set of values and 0.0 is returned for a &lt;= 3 value set.
*/
public double getKurtosis() {
* @return The kurtosis of a set of values. Double.NaN is returned for
* an empty set of values and 0.0 is returned for a &lt;= 3 value set.
*/
public double getKurtosis() {
if( n < 1) return Double.NaN;
if( n <= 3 ) return 0.0;
if( n < 1) return Double.NaN;
if( n <= 3 ) return 0.0;
double x1 = -6*Math.pow(sum,4);
double x2 = 12*((double)n)*Math.pow(sum,2)*sumsq;
double x3 = -3*((double)(n*(n-1)))*Math.pow(sumsq,2);
double x4 = -4*((double)(n*(n+1)))*sum*sumCube;
double x5 = Math.pow(((double)n),2)*((double)(n+1))*sumQuad;
double x1 = -6 * Math.pow(sum, 4);
double x2 = 12 * ((double) n) * Math.pow(sum, 2) * sumsq;
double x3 = -3 * ((double) (n * (n - 1))) * Math.pow(sumsq,2);
double x4 = -4 * ((double) (n * (n + 1))) * sum * sumCube;
double x5 = Math.pow(((double) n),2) * ((double) (n+1)) * sumQuad;
return (x1 + x2 + x3 + x4 + x5) /
( (double)(n*(n-1)*(n-2)*(n-3)) );
}
return (x1 + x2 + x3 + x4 + x5) /
( (double) (n * (n - 1) * (n - 2) * (n - 3)) );
}
/**
* Called in "addValue" to insert a new value into the statistic.
* @param v The value to be added.
*/
private void insertValue(double v) {
* @param v The value to be added.
*/
private void insertValue(double v) {
// The default value of product is NaN, if you
// try to retrieve the product for a univariate with
@ -250,34 +250,36 @@ public class UnivariateImpl implements Univariate, Serializable {
// Remove the influence of the discarded
sum -= discarded;
sumsq -= discarded * discarded;
sumCube -= Math.pow(discarded,3);
sumQuad -= Math.pow(discarded,4);
sumCube -= Math.pow(discarded, 3);
sumQuad -= Math.pow(discarded, 4);
if(discarded == min) {
min = doubleArray.getMin();
} else {
if(discarded == max){
} else if(discarded == max){
max = doubleArray.getMax();
}
}
if(product != 0.0){
// can safely remove discarded value
product *= v/discarded;
product *= v / discarded;
} else if(discarded == 0.0){
// need to recompute product
product = 1.0;
double[] elements = doubleArray.getElements();
for( int i = 0; i < elements.length; i++ ) {
product *= elements[i];
product *= elements[i];
}
} // else product = 0 and will still be 0 after discard
} else {
doubleArray.addElement( v );
n += 1.0;
if (v < min) min = v;
if (v > max) max = v;
if (v < min) {
min = v;
}
if (v > max) {
max = v;
}
product *= v;
}
} else {
@ -285,15 +287,19 @@ public class UnivariateImpl implements Univariate, Serializable {
// worry about storing any values. We don't need to discard the
// influence of any single item.
n += 1.0;
if (v < min) min = v;
if (v > max) max = v;
if (v < min) {
min = v;
}
if (v > max) {
max = v;
}
product *= v;
}
sum += v;
sumsq += v*v;
sumCube += Math.pow(v,3);
sumQuad += Math.pow(v,4);
sum += v;
sumsq += v * v;
sumCube += Math.pow(v,3);
sumQuad += Math.pow(v,4);
}
/** Getter for property max.
@ -339,19 +345,19 @@ public class UnivariateImpl implements Univariate, Serializable {
return sumsq;
}
/** Getter for property sumCube.
* @return Value of property sumCube.
*/
public double getSumCube() {
return sumCube;
}
/** Getter for property sumCube.
* @return Value of property sumCube.
*/
public double getSumCube() {
return sumCube;
}
/** Getter for property sumQuad.
* @return Value of property sumQuad.
*/
public double getSumQuad() {
return sumQuad;
}
/** Getter for property sumQuad.
* @return Value of property sumQuad.
*/
public double getSumQuad() {
return sumQuad;
}
/**
* Generates a text report displaying
@ -367,8 +373,8 @@ public class UnivariateImpl implements Univariate, Serializable {
outBuffer.append("max: " + max + "\n");
outBuffer.append("mean: " + getMean() + "\n");
outBuffer.append("std dev: " + getStandardDeviation() + "\n");
outBuffer.append("skewness: " + getSkewness() + "\n");
outBuffer.append("kurtosis: " + getKurtosis() + "\n");
outBuffer.append("skewness: " + getSkewness() + "\n");
outBuffer.append("kurtosis: " + getKurtosis() + "\n");
return outBuffer.toString();
}

88
src/test/org/apache/commons/math/distribution/TDistributionTest.java
View File

@ -85,45 +85,93 @@ public class TDistributionTest extends TestCase {
super.tearDown();
}
public void testLowerTailProbability(){
testProbability(-5.893, .001);
testProbability(-3.365, .010);
testProbability(-2.571, .025);
testProbability(-2.015, .050);
testProbability(-1.476, .100);
}
public void testUpperTailProbability(){
testProbability(5.893, .999);
testProbability(3.365, .990);
testProbability(2.571, .975);
testProbability(2.015, .950);
testProbability(1.476, .900);
}
public void testLowerTailValues(){
public void testInverseCummulativeProbability001() {
testValue(-5.893, .001);
}
public void testInverseCumulativeProbability010() {
testValue(-3.365, .010);
}
public void testInverseCumulativeProbability025() {
testValue(-2.571, .025);
}
public void testInverseCumulativeProbability050() {
testValue(-2.015, .050);
}
public void testInverseCumulativeProbability100() {
testValue(-1.476, .100);
}
public void testUpperTailValues(){
public void testInverseCummulativeProbability999() {
testValue(5.893, .999);
}
public void testInverseCumulativeProbability990() {
testValue(3.365, .990);
}
public void testInverseCumulativeProbability975() {
testValue(2.571, .975);
}
public void testInverseCumulativeProbability950() {
testValue(2.015, .950);
}
public void testInverseCumulativeProbability900() {
testValue(1.476, .900);
}
public void testCummulativeProbability001() {
testProbability(-5.893, .001);
}
public void testCumulativeProbability010() {
testProbability(-3.365, .010);
}
public void testCumulativeProbability025() {
testProbability(-2.571, .025);
}
public void testCumulativeProbability050() {
testProbability(-2.015, .050);
}
public void testCumulativeProbability100() {
testProbability(-1.476, .100);
}
public void testCummulativeProbability999() {
testProbability(5.893, .999);
}
public void testCumulativeProbability990() {
testProbability(3.365, .990);
}
public void testCumulativeProbability975() {
testProbability(2.571, .975);
}
public void testCumulativeProbability950() {
testProbability(2.015, .950);
}
public void testCumulativeProbability900() {
testProbability(1.476, .900);
}
private void testProbability(double x, double expected){
double actual = t.cummulativeProbability(x);
assertEquals("probability for " + x, expected, actual, 10e-4);
assertEquals(expected, actual, 10e-4);
}
private void testValue(double expected, double p){
double actual = t.inverseCummulativeProbability(p);
assertEquals("value for " + p, expected, actual, 10e-4);
assertEquals(expected, actual, 10e-4);
}
}

155
src/test/org/apache/commons/math/special/BetaTest.java Normal file
View File

@ -0,0 +1,155 @@
/* ====================================================================
* The Apache Software License, Version 1.1
*
* Copyright (c) 2003 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 acknowlegement:
* "This product includes software developed by the
* Apache Software Foundation (http://www.apache.org/)."
* Alternately, this acknowlegement may appear in the software itself,
* if and wherever such third-party acknowlegements 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 names 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.special;
import junit.framework.TestCase;
/**
* @author Brent Worden
*/
public class BetaTest extends TestCase {
/**
* Constructor for BetaTest.
* @param name
*/
public BetaTest(String name) {
super(name);
}
private void testRegularizedBeta(double expected, double x, double a, double b) {
double actual = Beta.regularizedBeta(x, a, b);
if (Double.isNaN(expected)) {
assertTrue(Double.isNaN(actual));
} else {
assertEquals(expected, actual, 10e-5);
}
}
private void testLogBeta(double expected, double a, double b) {
double actual = Beta.logBeta(a, b);
if (Double.isNaN(expected)) {
assertTrue(Double.isNaN(actual));
} else {
assertEquals(expected, actual, 10e-5);
}
}
public void testRegularizedBetaNanPositivePositive() {
testRegularizedBeta(Double.NaN, Double.NaN, 1.0, 1.0);
}
public void testRegularizedBetaPositiveNanPositive() {
testRegularizedBeta(Double.NaN, 0.5, Double.NaN, 1.0);
}
public void testRegularizedBetaPositivePositiveNan() {
testRegularizedBeta(Double.NaN, 0.5, 1.0, Double.NaN);
}
public void testRegularizedBetaNegativePositivePositive() {
testRegularizedBeta(Double.NaN, -0.5, 1.0, 2.0);
}
public void testRegularizedBetaPositiveNegativePositive() {
testRegularizedBeta(Double.NaN, 0.5, -1.0, 2.0);
}
public void testRegularizedBetaPositivePositiveNegative() {
testRegularizedBeta(Double.NaN, 0.5, 1.0, -2.0);
}
public void testRegularizedBetaZeroPositivePositive() {
testRegularizedBeta(0.0, 0.0, 1.0, 2.0);
}
public void testRegularizedBetaPositiveZeroPositive() {
testRegularizedBeta(Double.NaN, 0.5, 0.0, 2.0);
}
public void testRegularizedBetaPositivePositiveZero() {
testRegularizedBeta(Double.NaN, 0.5, 1.0, 0.0);
}
public void testRegularizedBetaPositivePositivePositive() {
testRegularizedBeta(0.75, 0.5, 1.0, 2.0);
}
public void testLogBetaNanPositive() {
testLogBeta(Double.NaN, Double.NaN, 2.0);
}
public void testLogBetaPositiveNan() {
testLogBeta(Double.NaN, 1.0, Double.NaN);
}
public void testLogBetaNegativePositive() {
testLogBeta(Double.NaN, -1.0, 2.0);
}
public void testLogBetaPositiveNegative() {
testLogBeta(Double.NaN, 1.0, -2.0);
}
public void testLogBetaZeroPositive() {
testLogBeta(Double.NaN, 0.0, 2.0);
}
public void testLogBetaPositiveZero() {
testLogBeta(Double.NaN, 1.0, 0.0);
}
public void testLogBetaPositivePositive() {
testLogBeta(-0.693147, 1.0, 2.0);
}
}