Apache/commons-math
1
0
Fork 0
mirror of https://github.com/apache/commons-math.git synced 2025-02-06 10:09:26 +00:00
Code Issues Packages Projects Releases Wiki Activity

- Merged ExponentialDistribution and ExponentialDistributionImpl (MATH-711).

Browse Source 
- Merged FDistribution and FDistributionImpl (MATH-711).
- Merged GammaDistribution and GammaDistributionImpl (MATH-711).

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1206399 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Sebastien Brisard 2011-11-26 06:17:49 +00:00
parent a050013f44
commit 32354a1039
13 changed files with 977 additions and 1195 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

428
src/main/java/org/apache/commons/math/distribution/ChiSquaredDistribution.java
View File

@ -1,214 +1,214 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
/**
* Implementation of the chi-squared distribution.
*
* @see <a href="http://en.wikipedia.org/wiki/Chi-squared_distribution">Chi-squared distribution (Wikipedia)</a>
* @see <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">Chi-squared Distribution (MathWorld)</a>
* @version $Id$
*/
public class ChiSquaredDistribution
extends AbstractContinuousDistribution
implements Serializable {
/**
* Default inverse cumulative probability accuracy
* @since 2.1
*/
public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
/** Serializable version identifier */
private static final long serialVersionUID = -8352658048349159782L;
/** Internal Gamma distribution. */
private final GammaDistribution gamma;
/** Inverse cumulative probability accuracy */
private final double solverAbsoluteAccuracy;
/**
* Create a Chi-Squared distribution with the given degrees of freedom.
*
* @param degreesOfFreedom Degrees of freedom.
*/
public ChiSquaredDistribution(double degreesOfFreedom) {
this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Create a Chi-Squared distribution with the given degrees of freedom and
* inverse cumulative probability accuracy.
*
* @param degreesOfFreedom Degrees of freedom.
* @param inverseCumAccuracy the maximum absolute error in inverse
* cumulative probability estimates (defaults to
* {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
* @since 2.1
*/
public ChiSquaredDistribution(double degreesOfFreedom,
double inverseCumAccuracy) {
gamma = new GammaDistributionImpl(degreesOfFreedom / 2, 2);
solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
* Access the number of degrees of freedom.
*
* @return the degrees of freedom.
*/
public double getDegreesOfFreedom() {
return gamma.getAlpha() * 2.0;
}
/** {@inheritDoc} */
public double density(double x) {
return gamma.density(x);
}
/** {@inheritDoc} */
public double cumulativeProbability(double x) {
return gamma.cumulativeProbability(x);
}
/**
* {@inheritDoc}
*
* Returns {@code 0} when {@code p == 0} and
* {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
*/
@Override
public double inverseCumulativeProbability(final double p) {
if (p == 0) {
return 0d;
}
if (p == 1) {
return Double.POSITIVE_INFINITY;
}
return super.inverseCumulativeProbability(p);
}
/** {@inheritDoc} */
@Override
protected double getDomainLowerBound(double p) {
return Double.MIN_VALUE * gamma.getBeta();
}
/** {@inheritDoc} */
@Override
protected double getDomainUpperBound(double p) {
// NOTE: chi squared is skewed to the left
// NOTE: therefore, P(X < &mu;) > .5
double ret;
if (p < .5) {
// use mean
ret = getDegreesOfFreedom();
} else {
// use max
ret = Double.MAX_VALUE;
}
return ret;
}
/** {@inheritDoc} */
@Override
protected double getInitialDomain(double p) {
// NOTE: chi squared is skewed to the left
// NOTE: therefore, P(X < &mu;) > 0.5
double ret;
if (p < 0.5) {
// use 1/2 mean
ret = getDegreesOfFreedom() * 0.5;
} else {
// use mean
ret = getDegreesOfFreedom();
}
return ret;
}
/** {@inheritDoc} */
@Override
protected double getSolverAbsoluteAccuracy() {
return solverAbsoluteAccuracy;
}
/**
* {@inheritDoc}
*
* The lower bound of the support is always 0 no matter the
* degrees of freedom.
*
* @return lower bound of the support (always 0)
*/
@Override
public double getSupportLowerBound() {
return 0;
}
/**
* {@inheritDoc}
*
* The upper bound of the support is always positive infinity no matter the
* degrees of freedom.
*
* @return upper bound of the support (always Double.POSITIVE_INFINITY)
*/
@Override
public double getSupportUpperBound() {
return Double.POSITIVE_INFINITY;
}
/**
* {@inheritDoc}
*
* For {@code k} degrees of freedom, the mean is {@code k}.
*/
@Override
protected double calculateNumericalMean() {
return getDegreesOfFreedom();
}
/**
* {@inheritDoc}
*
* For {@code k} degrees of freedom, the variance is {@code 2 * k}.
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalVariance() {
return 2*getDegreesOfFreedom();
}
/** {@inheritDoc} */
@Override
public boolean isSupportLowerBoundInclusive() {
return true;
}
/** {@inheritDoc} */
@Override
public boolean isSupportUpperBoundInclusive() {
return false;
}
}
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
/**
* Implementation of the chi-squared distribution.
*
* @see <a href="http://en.wikipedia.org/wiki/Chi-squared_distribution">Chi-squared distribution (Wikipedia)</a>
* @see <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">Chi-squared Distribution (MathWorld)</a>
* @version $Id: ChiSquaredDistribution.java 1206060 2011-11-25 05:16:56Z celestin $
*/
public class ChiSquaredDistribution
extends AbstractContinuousDistribution
implements Serializable {
/**
* Default inverse cumulative probability accuracy
* @since 2.1
*/
public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
/** Serializable version identifier */
private static final long serialVersionUID = -8352658048349159782L;
/** Internal Gamma distribution. */
private final GammaDistribution gamma;
/** Inverse cumulative probability accuracy */
private final double solverAbsoluteAccuracy;
/**
* Create a Chi-Squared distribution with the given degrees of freedom.
*
* @param degreesOfFreedom Degrees of freedom.
*/
public ChiSquaredDistribution(double degreesOfFreedom) {
this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Create a Chi-Squared distribution with the given degrees of freedom and
* inverse cumulative probability accuracy.
*
* @param degreesOfFreedom Degrees of freedom.
* @param inverseCumAccuracy the maximum absolute error in inverse
* cumulative probability estimates (defaults to
* {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
* @since 2.1
*/
public ChiSquaredDistribution(double degreesOfFreedom,
double inverseCumAccuracy) {
gamma = new GammaDistribution(degreesOfFreedom / 2, 2);
solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
* Access the number of degrees of freedom.
*
* @return the degrees of freedom.
*/
public double getDegreesOfFreedom() {
return gamma.getAlpha() * 2.0;
}
/** {@inheritDoc} */
public double density(double x) {
return gamma.density(x);
}
/** {@inheritDoc} */
public double cumulativeProbability(double x) {
return gamma.cumulativeProbability(x);
}
/**
* {@inheritDoc}
*
* Returns {@code 0} when {@code p == 0} and
* {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
*/
@Override
public double inverseCumulativeProbability(final double p) {
if (p == 0) {
return 0d;
}
if (p == 1) {
return Double.POSITIVE_INFINITY;
}
return super.inverseCumulativeProbability(p);
}
/** {@inheritDoc} */
@Override
protected double getDomainLowerBound(double p) {
return Double.MIN_VALUE * gamma.getBeta();
}
/** {@inheritDoc} */
@Override
protected double getDomainUpperBound(double p) {
// NOTE: chi squared is skewed to the left
// NOTE: therefore, P(X < &mu;) > .5
double ret;
if (p < .5) {
// use mean
ret = getDegreesOfFreedom();
} else {
// use max
ret = Double.MAX_VALUE;
}
return ret;
}
/** {@inheritDoc} */
@Override
protected double getInitialDomain(double p) {
// NOTE: chi squared is skewed to the left
// NOTE: therefore, P(X < &mu;) > 0.5
double ret;
if (p < 0.5) {
// use 1/2 mean
ret = getDegreesOfFreedom() * 0.5;
} else {
// use mean
ret = getDegreesOfFreedom();
}
return ret;
}
/** {@inheritDoc} */
@Override
protected double getSolverAbsoluteAccuracy() {
return solverAbsoluteAccuracy;
}
/**
* {@inheritDoc}
*
* The lower bound of the support is always 0 no matter the
* degrees of freedom.
*
* @return lower bound of the support (always 0)
*/
@Override
public double getSupportLowerBound() {
return 0;
}
/**
* {@inheritDoc}
*
* The upper bound of the support is always positive infinity no matter the
* degrees of freedom.
*
* @return upper bound of the support (always Double.POSITIVE_INFINITY)
*/
@Override
public double getSupportUpperBound() {
return Double.POSITIVE_INFINITY;
}
/**
* {@inheritDoc}
*
* For {@code k} degrees of freedom, the mean is {@code k}.
*/
@Override
protected double calculateNumericalMean() {
return getDegreesOfFreedom();
}
/**
* {@inheritDoc}
*
* For {@code k} degrees of freedom, the variance is {@code 2 * k}.
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalVariance() {
return 2*getDegreesOfFreedom();
}
/** {@inheritDoc} */
@Override
public boolean isSupportLowerBoundInclusive() {
return true;
}
/** {@inheritDoc} */
@Override
public boolean isSupportUpperBoundInclusive() {
return false;
}
}

232
src/main/java/org/apache/commons/math/distribution/ExponentialDistribution.java
View File

@ -16,24 +16,234 @@
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.FastMath;
/**
* The Exponential Distribution.
*
* <p>
* References:
* <ul>
* <li><a href="http://mathworld.wolfram.com/ExponentialDistribution.html">
* Exponential Distribution</a></li>
* </ul>
* </p>
* Implementation of the exponential distribution.
*
* @see <a href="http://en.wikipedia.org/wiki/Exponential_distribution">Exponential distribution (Wikipedia)</a>
* @see <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">Exponential distribution (MathWorld)</a>
* @version $Id$
*/
public interface ExponentialDistribution extends ContinuousDistribution {
public class ExponentialDistribution extends AbstractContinuousDistribution
implements Serializable {
/**
* Default inverse cumulative probability accuracy.
* @since 2.1
*/
public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
/** Serializable version identifier */
private static final long serialVersionUID = 2401296428283614780L;
/** The mean of this distribution. */
private final double mean;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
/**
* Create a exponential distribution with the given mean.
* @param mean mean of this distribution.
*/
public ExponentialDistribution(double mean) {
this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Create a exponential distribution with the given mean.
*
* @param mean Mean of this distribution.
* @param inverseCumAccuracy Maximum absolute error in inverse
* cumulative probability estimates (defaults to
* {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
* @throws NotStrictlyPositiveException if {@code mean <= 0}.
* @since 2.1
*/
public ExponentialDistribution(double mean, double inverseCumAccuracy)
throws NotStrictlyPositiveException{
if (mean <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
}
this.mean = mean;
solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
* Access the mean.
*
* @return the mean.
*/
double getMean();
public double getMean() {
return mean;
}
/** {@inheritDoc} */
public double density(double x) {
if (x < 0) {
return 0;
}
return FastMath.exp(-x / mean) / mean;
}
/**
* {@inheritDoc}
*
* The implementation of this method is based on:
* <ul>
* <li>
* <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">
* Exponential Distribution</a>, equation (1).</li>
* </ul>
*/
public double cumulativeProbability(double x) {
double ret;
if (x <= 0.0) {
ret = 0.0;
} else {
ret = 1.0 - FastMath.exp(-x / mean);
}
return ret;
}
/**
* {@inheritDoc}
*
* Returns {@code 0} when {@code p= = 0} and
* {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
*/
@Override
public double inverseCumulativeProbability(double p) throws OutOfRangeException {
double ret;
if (p < 0.0 || p > 1.0) {
throw new OutOfRangeException(p, 0.0, 1.0);
} else if (p == 1.0) {
ret = Double.POSITIVE_INFINITY;
} else {
ret = -mean * FastMath.log(1.0 - p);
}
return ret;
}
/**
* {@inheritDoc}
*
* <p><strong>Algorithm Description</strong>: this implementation uses the
* <a href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html">
* Inversion Method</a> to generate exponentially distributed random values
* from uniform deviates.</p>
*
* @return a random value.
* @since 2.2
*/
@Override
public double sample() {
return randomData.nextExponential(mean);
}
/** {@inheritDoc} */
@Override
protected double getDomainLowerBound(double p) {
return 0;
}
/** {@inheritDoc} */
@Override
protected double getDomainUpperBound(double p) {
// NOTE: exponential is skewed to the left
// NOTE: therefore, P(X < &mu;) > .5
if (p < 0.5) {
// use mean
return mean;
} else {
// use max
return Double.MAX_VALUE;
}
}
/** {@inheritDoc} */
@Override
protected double getInitialDomain(double p) {
// TODO: try to improve on this estimate
// TODO: what should really happen here is not derive from
// AbstractContinuousDistribution
// TODO: because the inverse cumulative distribution is simple.
// Exponential is skewed to the left, therefore, P(X < &mu;) > .5
if (p < 0.5) {
// use 1/2 mean
return mean * 0.5;
} else {
// use mean
return mean;
}
}
/** {@inheritDoc} */
@Override
protected double getSolverAbsoluteAccuracy() {
return solverAbsoluteAccuracy;
}
/**
* {@inheritDoc}
*
* The lower bound of the support is always 0 no matter the mean parameter.
*
* @return lower bound of the support (always 0)
*/
@Override
public double getSupportLowerBound() {
return 0;
}
/**
* {@inheritDoc}
*
* The upper bound of the support is always positive infinity
* no matter the mean parameter.
*
* @return upper bound of the support (always Double.POSITIVE_INFINITY)
*/
@Override
public double getSupportUpperBound() {
return Double.POSITIVE_INFINITY;
}
/**
* {@inheritDoc}
*
* For mean parameter {@code k}, the mean is {@code k}.
*/
@Override
protected double calculateNumericalMean() {
return getMean();
}
/**
* {@inheritDoc}
*
* For mean parameter {@code k}, the variance is {@code k^2}.
*/
@Override
protected double calculateNumericalVariance() {
final double m = getMean();
return m * m;
}
/** {@inheritDoc} */
@Override
public boolean isSupportLowerBoundInclusive() {
return true;
}
/** {@inheritDoc} */
@Override
public boolean isSupportUpperBoundInclusive() {
return false;
}
}

279
src/main/java/org/apache/commons/math/distribution/ExponentialDistributionImpl.java
View File

@ -1,279 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.FastMath;
/**
* The default implementation of {@link ExponentialDistribution}.
*
* @version $Id$
*/
public class ExponentialDistributionImpl extends AbstractContinuousDistribution
implements ExponentialDistribution, Serializable {
/**
* Default inverse cumulative probability accuracy.
* @since 2.1
*/
public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
/** Serializable version identifier */
private static final long serialVersionUID = 2401296428283614780L;
/** The mean of this distribution. */
private final double mean;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
/**
* Create a exponential distribution with the given mean.
* @param mean mean of this distribution.
*/
public ExponentialDistributionImpl(double mean) {
this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Create a exponential distribution with the given mean.
*
* @param mean Mean of this distribution.
* @param inverseCumAccuracy Maximum absolute error in inverse
* cumulative probability estimates (defaults to
* {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
* @throws NotStrictlyPositiveException if {@code mean <= 0}.
* @since 2.1
*/
public ExponentialDistributionImpl(double mean, double inverseCumAccuracy) {
if (mean <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
}
this.mean = mean;
solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
* {@inheritDoc}
*/
public double getMean() {
return mean;
}
/**
* {@inheritDoc}
*/
public double density(double x) {
if (x < 0) {
return 0;
}
return FastMath.exp(-x / mean) / mean;
}
/**
* {@inheritDoc}
*
* The implementation of this method is based on:
* <ul>
* <li>
* <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">
* Exponential Distribution</a>, equation (1).</li>
* </ul>
*/
public double cumulativeProbability(double x) {
double ret;
if (x <= 0.0) {
ret = 0.0;
} else {
ret = 1.0 - FastMath.exp(-x / mean);
}
return ret;
}
/**
* {@inheritDoc}
*
* It will return {@code 0} when {@code p = 0} and
* {@code Double.POSITIVE_INFINITY} when {@code p = 1}.
*/
@Override
public double inverseCumulativeProbability(double p) throws OutOfRangeException {
double ret;
if (p < 0.0 || p > 1.0) {
throw new OutOfRangeException(p, 0.0, 1.0);
} else if (p == 1.0) {
ret = Double.POSITIVE_INFINITY;
} else {
ret = -mean * FastMath.log(1.0 - p);
}
return ret;
}
/**
* Generates a random value sampled from this distribution.
*
* <p><strong>Algorithm Description</strong>: Uses the <a
* href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html"> Inversion
* Method</a> to generate exponentially distributed random values from
* uniform deviates.</p>
*
* @return a random value.
* @since 2.2
*/
@Override
public double sample() {
return randomData.nextExponential(mean);
}
/**
* Access the domain value lower bound, based on {@code p}, used to
* bracket a CDF root.
*
* @param p Desired probability for the critical value.
* @return the domain value lower bound, i.e. {@code P(X < 'lower bound') < p}.
*/
@Override
protected double getDomainLowerBound(double p) {
return 0;
}
/**
* Access the domain value upper bound, based on {@code p}, used to
* bracket a CDF root.
*
* @param p Desired probability for the critical value.
* @return the domain value upper bound, i.e. {@code P(X < 'upper bound') > p}.
*/
@Override
protected double getDomainUpperBound(double p) {
// NOTE: exponential is skewed to the left
// NOTE: therefore, P(X < &mu;) > .5
if (p < 0.5) {
// use mean
return mean;
} else {
// use max
return Double.MAX_VALUE;
}
}
/**
* Access the initial domain value, based on {@code p}, used to
* bracket a CDF root.
*
* @param p Desired probability for the critical value.
* @return the initial domain value.
*/
@Override
protected double getInitialDomain(double p) {
// TODO: try to improve on this estimate
// TODO: what should really happen here is not derive from AbstractContinuousDistribution
// TODO: because the inverse cumulative distribution is simple.
// Exponential is skewed to the left, therefore, P(X < &mu;) > .5
if (p < 0.5) {
// use 1/2 mean
return mean * 0.5;
} else {
// use mean
return mean;
}
}
/**
* Return the absolute accuracy setting of the solver used to estimate
* inverse cumulative probabilities.
*
* @return the solver absolute accuracy.
* @since 2.1
*/
@Override
protected double getSolverAbsoluteAccuracy() {
return solverAbsoluteAccuracy;
}
/**
* {@inheritDoc}
*
* The lower bound of the support is always 0 no matter the mean parameter.
*
* @return lower bound of the support (always 0)
*/
@Override
public double getSupportLowerBound() {
return 0;
}
/**
* {@inheritDoc}
*
* The upper bound of the support is always positive infinity
* no matter the mean parameter.
*
* @return upper bound of the support (always Double.POSITIVE_INFINITY)
*/
@Override
public double getSupportUpperBound() {
return Double.POSITIVE_INFINITY;
}
/**
* {@inheritDoc}
*
* For mean parameter <code>k</code>, the mean is
* <code>k</code>
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalMean() {
return getMean();
}
/**
* {@inheritDoc}
*
* For mean parameter <code>k</code>, the variance is
* <code>k^2</code>
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalVariance() {
final double m = getMean();
return m * m;
}
/**
* {@inheritDoc}
*/
@Override
public boolean isSupportLowerBoundInclusive() {
return true;
}
/**
* {@inheritDoc}
*/
@Override
public boolean isSupportUpperBoundInclusive() {
return false;
}
}

268
src/main/java/org/apache/commons/math/distribution/FDistribution.java
View File

@ -14,33 +14,277 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.special.Beta;
import org.apache.commons.math.util.FastMath;
/**
* F-Distribution.
*
* <p>
* References:
* <ul>
* <li><a href="http://mathworld.wolfram.com/F-Distribution.html">
* F-Distribution</a></li>
* </ul>
* </p>
* Implementation of the F-distribution.
*
* @see <a href="http://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a>
* @see <a href="http://mathworld.wolfram.com/F-Distribution.html">F-distribution (MathWorld)</a>
* @version $Id$
*/
public interface FDistribution extends ContinuousDistribution {
public class FDistribution
extends AbstractContinuousDistribution
implements Serializable {
/**
* Default inverse cumulative probability accuracy.
* @since 2.1
*/
public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
/** Serializable version identifier. */
private static final long serialVersionUID = -8516354193418641566L;
/** The numerator degrees of freedom. */
private final double numeratorDegreesOfFreedom;
/** The numerator degrees of freedom. */
private final double denominatorDegreesOfFreedom;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
/**
* Create a F distribution using the given degrees of freedom.
* @param numeratorDegreesOfFreedom Numerator degrees of freedom.
* @param denominatorDegreesOfFreedom Denominator degrees of freedom.
* @throws NotStrictlyPositiveException if
* {@code numeratorDegreesOfFreedom <= 0} or
* {@code denominatorDegreesOfFreedom <= 0}.
*/
public FDistribution(double numeratorDegreesOfFreedom,
double denominatorDegreesOfFreedom)
throws NotStrictlyPositiveException {
this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom,
DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Create an F distribution using the given degrees of freedom
* and inverse cumulative probability accuracy.
* @param numeratorDegreesOfFreedom Numerator degrees of freedom.
* @param denominatorDegreesOfFreedom Denominator degrees of freedom.
* @param inverseCumAccuracy the maximum absolute error in inverse
* cumulative probability estimates.
* (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
* @throws NotStrictlyPositiveException if
* {@code numeratorDegreesOfFreedom <= 0} or
* {@code denominatorDegreesOfFreedom <= 0}.
* @since 2.1
*/
public FDistribution(double numeratorDegreesOfFreedom,
double denominatorDegreesOfFreedom,
double inverseCumAccuracy)
throws NotStrictlyPositiveException {
if (numeratorDegreesOfFreedom <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
numeratorDegreesOfFreedom);
}
if (denominatorDegreesOfFreedom <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
denominatorDegreesOfFreedom);
}
this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
* {@inheritDoc}
*
* @since 2.1
*/
public double density(double x) {
final double nhalf = numeratorDegreesOfFreedom / 2;
final double mhalf = denominatorDegreesOfFreedom / 2;
final double logx = FastMath.log(x);
final double logn = FastMath.log(numeratorDegreesOfFreedom);
final double logm = FastMath.log(denominatorDegreesOfFreedom);
final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
denominatorDegreesOfFreedom);
return FastMath.exp(nhalf * logn + nhalf * logx - logx +
mhalf * logm - nhalf * lognxm - mhalf * lognxm -
Beta.logBeta(nhalf, mhalf));
}
/**
* {@inheritDoc}
*
* The implementation of this method is based on
* <ul>
* <li>
* <a href="http://mathworld.wolfram.com/F-Distribution.html">
* F-Distribution</a>, equation (4).
* </li>
* </ul>
*/
public double cumulativeProbability(double x) {
double ret;
if (x <= 0) {
ret = 0;
} else {
double n = numeratorDegreesOfFreedom;
double m = denominatorDegreesOfFreedom;
ret = Beta.regularizedBeta((n * x) / (m + n * x),
0.5 * n,
0.5 * m);
}
return ret;
}
/**
* {@inheritDoc}
*
* Returns {@code 0} when {@code p == 0} and
* {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
*/
@Override
public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
if (p == 0) {
return 0;
}
if (p == 1) {
return Double.POSITIVE_INFINITY;
}
return super.inverseCumulativeProbability(p);
}
/** {@inheritDoc} */
@Override
protected double getDomainLowerBound(double p) {
return 0;
}
/** {@inheritDoc} */
@Override
protected double getDomainUpperBound(double p) {
return Double.MAX_VALUE;
}
/** {@inheritDoc} */
@Override
protected double getInitialDomain(double p) {
double ret = 1;
double d = denominatorDegreesOfFreedom;
if (d > 2) {
// use mean
ret = d / (d - 2);
}
return ret;
}
/**
* Access the numerator degrees of freedom.
*
* @return the numerator degrees of freedom.
*/
double getNumeratorDegreesOfFreedom();
public double getNumeratorDegreesOfFreedom() {
return numeratorDegreesOfFreedom;
}
/**
* Access the denominator degrees of freedom.
*
* @return the denominator degrees of freedom.
*/
double getDenominatorDegreesOfFreedom();
public double getDenominatorDegreesOfFreedom() {
return denominatorDegreesOfFreedom;
}
/** {@inheritDoc} */
@Override
protected double getSolverAbsoluteAccuracy() {
return solverAbsoluteAccuracy;
}
/**
* {@inheritDoc}
*
* The lower bound of the support is always 0 no matter the parameters.
*
* @return lower bound of the support (always 0)
*/
@Override
public double getSupportLowerBound() {
return 0;
}
/**
* {@inheritDoc}
*
* The upper bound of the support is always positive infinity
* no matter the parameters.
*
* @return upper bound of the support (always Double.POSITIVE_INFINITY)
*/
@Override
public double getSupportUpperBound() {
return Double.POSITIVE_INFINITY;
}
/**
* {@inheritDoc}
*
* For denominator degrees of freedom parameter {@code b}, the mean is
* <ul>
* <li>if {@code b > 2} then {@code b / (b - 2)},</li>
* <li>else undefined ({@code Double.NaN}).
* </ul>
*/
@Override
protected double calculateNumericalMean() {
final double denominatorDF = getDenominatorDegreesOfFreedom();
if (denominatorDF > 2) {
return denominatorDF / (denominatorDF - 2);
}
return Double.NaN;
}
/**
* {@inheritDoc}
*
* For numerator degrees of freedom parameter {@code a} and denominator
* degrees of freedom parameter {@code b}, the variance is
* <ul>
* <li>
* if {@code b > 4} then
* {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]},
* </li>
* <li>else undefined ({@code Double.NaN}).
* </ul>
*/
@Override
protected double calculateNumericalVariance() {
final double denominatorDF = getDenominatorDegreesOfFreedom();
if (denominatorDF > 4) {
final double numeratorDF = getNumeratorDegreesOfFreedom();
final double denomDFMinusTwo = denominatorDF - 2;
return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
}
return Double.NaN;
}
/** {@inheritDoc} */
@Override
public boolean isSupportLowerBoundInclusive() {
return true;
}
/** {@inheritDoc} */
@Override
public boolean isSupportUpperBoundInclusive() {
return false;
}
}

318
src/main/java/org/apache/commons/math/distribution/FDistributionImpl.java
View File

@ -1,318 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.special.Beta;
import org.apache.commons.math.util.FastMath;
/**
* Default implementation of
* {@link org.apache.commons.math.distribution.FDistribution}.
*
* @version $Id$
*/
public class FDistributionImpl
extends AbstractContinuousDistribution
implements FDistribution, Serializable {
/**
* Default inverse cumulative probability accuracy.
* @since 2.1
*/
public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
/** Serializable version identifier. */
private static final long serialVersionUID = -8516354193418641566L;
/** The numerator degrees of freedom. */
private final double numeratorDegreesOfFreedom;
/** The numerator degrees of freedom. */
private final double denominatorDegreesOfFreedom;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
/**
* Create a F distribution using the given degrees of freedom.
* @param numeratorDegreesOfFreedom Numerator degrees of freedom.
* @param denominatorDegreesOfFreedom Denominator degrees of freedom.
* @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0}
* or {@code denominatorDegreesOfFreedom <= 0}.
*/
public FDistributionImpl(double numeratorDegreesOfFreedom,
double denominatorDegreesOfFreedom) {
this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom,
DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Create an F distribution using the given degrees of freedom
* and inverse cumulative probability accuracy.
* @param numeratorDegreesOfFreedom Numerator degrees of freedom.
* @param denominatorDegreesOfFreedom Denominator degrees of freedom.
* @param inverseCumAccuracy the maximum absolute error in inverse
* cumulative probability estimates.
* (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
* @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0}
* or {@code denominatorDegreesOfFreedom <= 0}.
* @since 2.1
*/
public FDistributionImpl(double numeratorDegreesOfFreedom,
double denominatorDegreesOfFreedom,
double inverseCumAccuracy) {
if (numeratorDegreesOfFreedom <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
numeratorDegreesOfFreedom);
}
if (denominatorDegreesOfFreedom <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
denominatorDegreesOfFreedom);
}
this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
* {@inheritDoc}
*
* @since 2.1
*/
public double density(double x) {
final double nhalf = numeratorDegreesOfFreedom / 2;
final double mhalf = denominatorDegreesOfFreedom / 2;
final double logx = FastMath.log(x);
final double logn = FastMath.log(numeratorDegreesOfFreedom);
final double logm = FastMath.log(denominatorDegreesOfFreedom);
final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
denominatorDegreesOfFreedom);
return FastMath.exp(nhalf * logn + nhalf * logx - logx +
mhalf * logm - nhalf * lognxm - mhalf * lognxm -
Beta.logBeta(nhalf, mhalf));
}
/**
* {@inheritDoc}
*
* The implementation of this method is based on
* <ul>
* <li>
* <a href="http://mathworld.wolfram.com/F-Distribution.html">
* F-Distribution</a>, equation (4).
* </li>
* </ul>
*/
public double cumulativeProbability(double x) {
double ret;
if (x <= 0) {
ret = 0;
} else {
double n = numeratorDegreesOfFreedom;
double m = denominatorDegreesOfFreedom;
ret = Beta.regularizedBeta((n * x) / (m + n * x),
0.5 * n,
0.5 * m);
}
return ret;
}
/**
* {@inheritDoc}
*
* It will return {@code 0} when {@code p = 0} and
* {@code Double.POSITIVE_INFINITY} when {@code p = 1}.
*/
@Override
public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
if (p == 0) {
return 0;
}
if (p == 1) {
return Double.POSITIVE_INFINITY;
}
return super.inverseCumulativeProbability(p);
}
/**
* Access the domain value lower bound, based on {@code p}, used to
* bracket a CDF root. This method is used by
* {@link #inverseCumulativeProbability(double)} to find critical values.
*
* @param p Desired probability for the critical value.
* @return the domain value lower bound, i.e. {@code P(X < 'lower bound') < p}.
*/
@Override
protected double getDomainLowerBound(double p) {
return 0;
}
/**
* Access the domain value upper bound, based on {@code p}, used to
* bracket a CDF root. This method is used by
* {@link #inverseCumulativeProbability(double)} to find critical values.
*
* @param p Desired probability for the critical value.
* @return the domain value upper bound, i.e. {@code P(X < 'upper bound') > p}.
*/
@Override
protected double getDomainUpperBound(double p) {
return Double.MAX_VALUE;
}
/**
* Access the initial domain value, based on {@code p}, used to
* bracket a CDF root. This method is used by
* {@link #inverseCumulativeProbability(double)} to find critical values.
*
* @param p Desired probability for the critical value.
* @return the initial domain value.
*/
@Override
protected double getInitialDomain(double p) {
double ret = 1;
double d = denominatorDegreesOfFreedom;
if (d > 2) {
// use mean
ret = d / (d - 2);
}
return ret;
}
/**
* {@inheritDoc}
*/
public double getNumeratorDegreesOfFreedom() {
return numeratorDegreesOfFreedom;
}
/**
* {@inheritDoc}
*/
public double getDenominatorDegreesOfFreedom() {
return denominatorDegreesOfFreedom;
}
/**
* Return the absolute accuracy setting of the solver used to estimate
* inverse cumulative probabilities.
*
* @return the solver absolute accuracy
* @since 2.1
*/
@Override
protected double getSolverAbsoluteAccuracy() {
return solverAbsoluteAccuracy;
}
/**
* {@inheritDoc}
*
* The lower bound of the support is always 0 no matter the parameters.
*
* @return lower bound of the support (always 0)
*/
@Override
public double getSupportLowerBound() {
return 0;
}
/**
* {@inheritDoc}
*
* The upper bound of the support is always positive infinity
* no matter the parameters.
*
* @return upper bound of the support (always Double.POSITIVE_INFINITY)
*/
@Override
public double getSupportUpperBound() {
return Double.POSITIVE_INFINITY;
}
/**
* {@inheritDoc}
*
* For denominator degrees of freedom parameter <code>b</code>,
* the mean is
* <ul>
* <li>if <code>b &gt; 2</code> then <code>b / (b - 2)</code></li>
* <li>else <code>undefined</code>
* </ul>
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalMean() {
final double denominatorDF = getDenominatorDegreesOfFreedom();
if (denominatorDF > 2) {
return denominatorDF / (denominatorDF - 2);
}
return Double.NaN;
}
/**
* {@inheritDoc}
*
* For numerator degrees of freedom parameter <code>a</code>
* and denominator degrees of freedom parameter <code>b</code>,
* the variance is
* <ul>
* <li>
* if <code>b &gt; 4</code> then
* <code>[ 2 * b^2 * (a + b - 2) ] / [ a * (b - 2)^2 * (b - 4) ]</code>
* </li>
* <li>else <code>undefined</code>
* </ul>
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalVariance() {
final double denominatorDF = getDenominatorDegreesOfFreedom();
if (denominatorDF > 4) {
final double numeratorDF = getNumeratorDegreesOfFreedom();
final double denomDFMinusTwo = denominatorDF - 2;
return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
}
return Double.NaN;
}
/**
* {@inheritDoc}
*/
@Override
public boolean isSupportLowerBoundInclusive() {
return true;
}
/**
* {@inheritDoc}
*/
@Override
public boolean isSupportUpperBoundInclusive() {
return false;
}
}

259
src/main/java/org/apache/commons/math/distribution/GammaDistribution.java
View File

@ -16,31 +16,254 @@
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.special.Gamma;
import org.apache.commons.math.util.FastMath;
/**
* The Gamma Distribution.
*
* <p>
* References:
* <ul>
* <li><a href="http://mathworld.wolfram.com/GammaDistribution.html">
* Gamma Distribution</a></li>
* </ul>
* </p>
* Implementation of the Gamma distribution.
*
* @see <a href="http://en.wikipedia.org/wiki/Gamma_distribution">Gamma distribution (Wikipedia)</a>
* @see <a href="http://mathworld.wolfram.com/GammaDistribution.html">Gamma distribution (MathWorld)</a>
* @version $Id$
*/
public interface GammaDistribution extends ContinuousDistribution {
public class GammaDistribution extends AbstractContinuousDistribution
implements Serializable {
/**
* Access the alpha shape parameter.
*
* @return alpha.
* Default inverse cumulative probability accuracy.
* @since 2.1
*/
double getAlpha();
public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
/** Serializable version identifier. */
private static final long serialVersionUID = -3239549463135430361L;
/** The shape parameter. */
private final double alpha;
/** The scale parameter. */
private final double beta;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
/**
* Access the beta scale parameter.
*
* @return beta.
* Create a new gamma distribution with the given {@code alpha} and
* {@code beta} values.
* @param alpha the shape parameter.
* @param beta the scale parameter.
*/
double getBeta();
public GammaDistribution(double alpha, double beta) {
this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Create a new gamma distribution with the given {@code alpha} and
* {@code beta} values.
*
* @param alpha Shape parameter.
* @param beta Scale parameter.
* @param inverseCumAccuracy Maximum absolute error in inverse
* cumulative probability estimates (defaults to
* {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
* @throws NotStrictlyPositiveException if {@code alpha <= 0} or
* {@code beta <= 0}.
* @since 2.1
*/
public GammaDistribution(double alpha, double beta, double inverseCumAccuracy)
throws NotStrictlyPositiveException {
if (alpha <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.ALPHA, alpha);
}
if (beta <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.BETA, beta);
}
this.alpha = alpha;
this.beta = beta;
solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
* {@inheritDoc}
*
* The implementation of this method is based on:
* <ul>
* <li>
* <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
* Chi-Squared Distribution</a>, equation (9).
* </li>
* <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
* Belmont, CA: Duxbury Press.
* </li>
* </ul>
*/
public double cumulativeProbability(double x) {
double ret;
if (x <= 0) {
ret = 0;
} else {
ret = Gamma.regularizedGammaP(alpha, x / beta);
}
return ret;
}
/**
* {@inheritDoc}
*
* Returns {@code 0} when {@code p == 0} and
* {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
*/
@Override
public double inverseCumulativeProbability(final double p) {
if (p == 0) {
return 0;
}
if (p == 1) {
return Double.POSITIVE_INFINITY;
}
return super.inverseCumulativeProbability(p);
}
/**
* Access the {@code alpha} shape parameter.
*
* @return {@code alpha}.
*/
public double getAlpha() {
return alpha;
}
/**
* Access the {@code beta} scale parameter.
*
* @return {@code beta}.
*/
public double getBeta() {
return beta;
}
/** {@inheritDoc} */
public double density(double x) {
if (x < 0) {
return 0;
}
return FastMath.pow(x / beta, alpha - 1) / beta *
FastMath.exp(-x / beta) / FastMath.exp(Gamma.logGamma(alpha));
}
/** {@inheritDoc} */
@Override
protected double getDomainLowerBound(double p) {
// TODO: try to improve on this estimate
return Double.MIN_VALUE;
}
/** {@inheritDoc} */
@Override
protected double getDomainUpperBound(double p) {
// TODO: try to improve on this estimate
// NOTE: gamma is skewed to the left
// NOTE: therefore, P(X < &mu;) > .5
double ret;
if (p < 0.5) {
// use mean
ret = alpha * beta;
} else {
// use max value
ret = Double.MAX_VALUE;
}
return ret;
}
/** {@inheritDoc} */
@Override
protected double getInitialDomain(double p) {
// TODO: try to improve on this estimate
// Gamma is skewed to the left, therefore, P(X < &mu;) > .5
double ret;
if (p < 0.5) {
// use 1/2 mean
ret = alpha * beta * 0.5;
} else {
// use mean
ret = alpha * beta;
}
return ret;
}
/** {@inheritDoc} */
@Override
protected double getSolverAbsoluteAccuracy() {
return solverAbsoluteAccuracy;
}
/**
* {@inheritDoc}
*
* The lower bound of the support is always 0 no matter the parameters.
*
* @return lower bound of the support (always 0)
*/
@Override
public double getSupportLowerBound() {
return 0;
}
/**
* {@inheritDoc}
*
* The upper bound of the support is always positive infinity
* no matter the parameters.
*
* @return upper bound of the support (always Double.POSITIVE_INFINITY)
*/
@Override
public double getSupportUpperBound() {
return Double.POSITIVE_INFINITY;
}
/**
* {@inheritDoc}
*
* For shape parameter {@code alpha} and scale parameter {@code beta}, the
* mean is {@code alpha * beta}.
*/
@Override
protected double calculateNumericalMean() {
return getAlpha() * getBeta();
}
/**
* {@inheritDoc}
*
* For shape parameter {@code alpha} and scale parameter {@code beta}, the
* variance is {@code alpha * beta^2}.
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalVariance() {
final double b = getBeta();
return getAlpha() * b * b;
}
/** {@inheritDoc} */
@Override
public boolean isSupportLowerBoundInclusive() {
return true;
}
/** {@inheritDoc} */
@Override
public boolean isSupportUpperBoundInclusive() {
return false;
}
}

297
src/main/java/org/apache/commons/math/distribution/GammaDistributionImpl.java
View File

@ -1,297 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.special.Gamma;
import org.apache.commons.math.util.FastMath;
/**
* The default implementation of {@link GammaDistribution}.
*
* @version $Id$
*/
public class GammaDistributionImpl extends AbstractContinuousDistribution
implements GammaDistribution, Serializable {
/**
* Default inverse cumulative probability accuracy.
* @since 2.1
*/
public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
/** Serializable version identifier. */
private static final long serialVersionUID = -3239549463135430361L;
/** The shape parameter. */
private final double alpha;
/** The scale parameter. */
private final double beta;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
/**
* Create a new gamma distribution with the given alpha and beta values.
* @param alpha the shape parameter.
* @param beta the scale parameter.
*/
public GammaDistributionImpl(double alpha, double beta) {
this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Create a new gamma distribution with the given alpha and beta values.
*
* @param alpha Shape parameter.
* @param beta Scale parameter.
* @param inverseCumAccuracy Maximum absolute error in inverse
* cumulative probability estimates (defaults to
* {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
* @throws NotStrictlyPositiveException if {@code alpha <= 0} or
* {@code beta <= 0}.
* @since 2.1
*/
public GammaDistributionImpl(double alpha, double beta, double inverseCumAccuracy) {
if (alpha <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.ALPHA, alpha);
}
if (beta <= 0) {
throw new NotStrictlyPositiveException(LocalizedFormats.BETA, beta);
}
this.alpha = alpha;
this.beta = beta;
solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
* {@inheritDoc}
*
* The implementation of this method is based on:
* <ul>
* <li>
* <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
* Chi-Squared Distribution</a>, equation (9).
* </li>
* <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
* Belmont, CA: Duxbury Press.
* </li>
* </ul>
*/
public double cumulativeProbability(double x) {
double ret;
if (x <= 0) {
ret = 0;
} else {
ret = Gamma.regularizedGammaP(alpha, x / beta);
}
return ret;
}
/**
* {@inheritDoc}
*
* It will return {@code 0} when {@cod p = 0} and
* {@code Double.POSITIVE_INFINITY} when {@code p = 1}.
*/
@Override
public double inverseCumulativeProbability(final double p) {
if (p == 0) {
return 0;
}
if (p == 1) {
return Double.POSITIVE_INFINITY;
}
return super.inverseCumulativeProbability(p);
}
/**
* {@inheritDoc}
*/
public double getAlpha() {
return alpha;
}
/**
* {@inheritDoc}
*/
public double getBeta() {
return beta;
}
/**
* {@inheritDoc}
*/
public double density(double x) {
if (x < 0) {
return 0;
}
return FastMath.pow(x / beta, alpha - 1) / beta *
FastMath.exp(-x / beta) / FastMath.exp(Gamma.logGamma(alpha));
}
/**
* Access the domain value lower bound, based on {@code p}, used to
* bracket a CDF root. This method is used by
* {@link #inverseCumulativeProbability(double)} to find critical values.
*
* @param p Desired probability for the critical value.
* @return the domain value lower bound, i.e. {@code P(X < 'lower bound') < p}.
*/
@Override
protected double getDomainLowerBound(double p) {
// TODO: try to improve on this estimate
return Double.MIN_VALUE;
}
/**
* Access the domain value upper bound, based on {@code p}, used to
* bracket a CDF root. This method is used by
* {@link #inverseCumulativeProbability(double)} to find critical values.
*
* @param p Desired probability for the critical value.
* @return the domain value upper bound, i.e. {@code P(X < 'upper bound') > p}.
*/
@Override
protected double getDomainUpperBound(double p) {
// TODO: try to improve on this estimate
// NOTE: gamma is skewed to the left
// NOTE: therefore, P(X < &mu;) > .5
double ret;
if (p < 0.5) {
// use mean
ret = alpha * beta;
} else {
// use max value
ret = Double.MAX_VALUE;
}
return ret;
}
/**
* Access the initial domain value, based on {@code p}, used to
* bracket a CDF root. This method is used by
* {@link #inverseCumulativeProbability(double)} to find critical values.
*
* @param p Desired probability for the critical value.
* @return the initial domain value.
*/
@Override
protected double getInitialDomain(double p) {
// TODO: try to improve on this estimate
// Gamma is skewed to the left, therefore, P(X < &mu;) > .5
double ret;
if (p < 0.5) {
// use 1/2 mean
ret = alpha * beta * 0.5;
} else {
// use mean
ret = alpha * beta;
}
return ret;
}
/**
* Return the absolute accuracy setting of the solver used to estimate
* inverse cumulative probabilities.
*
* @return the solver absolute accuracy.
* @since 2.1
*/
@Override
protected double getSolverAbsoluteAccuracy() {
return solverAbsoluteAccuracy;
}
/**
* {@inheritDoc}
*
* The lower bound of the support is always 0 no matter the parameters.
*
* @return lower bound of the support (always 0)
*/
@Override
public double getSupportLowerBound() {
return 0;
}
/**
* {@inheritDoc}
*
* The upper bound of the support is always positive infinity
* no matter the parameters.
*
* @return upper bound of the support (always Double.POSITIVE_INFINITY)
*/
@Override
public double getSupportUpperBound() {
return Double.POSITIVE_INFINITY;
}
/**
* {@inheritDoc}
*
* For shape parameter <code>alpha</code> and scale
* parameter <code>beta</code>, the mean is
* <code>alpha * beta</code>
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalMean() {
return getAlpha() * getBeta();
}
/**
* {@inheritDoc}
*
* For shape parameter <code>alpha</code> and scale
* parameter <code>beta</code>, the variance is
* <code>alpha * beta^2</code>
*
* @return {@inheritDoc}
*/
@Override
protected double calculateNumericalVariance() {
final double b = getBeta();
return getAlpha() * b * b;
}
/**
* {@inheritDoc}
*/
@Override
public boolean isSupportLowerBoundInclusive() {
return true;
}
/**
* {@inheritDoc}
*/
@Override
public boolean isSupportUpperBoundInclusive() {
return false;
}
}

8
src/main/java/org/apache/commons/math/random/RandomDataImpl.java
View File

@ -29,7 +29,7 @@ import org.apache.commons.math.distribution.BinomialDistribution;
import org.apache.commons.math.distribution.CauchyDistribution;
import org.apache.commons.math.distribution.ChiSquaredDistribution;
import org.apache.commons.math.distribution.ContinuousDistribution;
import org.apache.commons.math.distribution.FDistributionImpl;
import org.apache.commons.math.distribution.FDistribution;
import org.apache.commons.math.distribution.HypergeometricDistributionImpl;
import org.apache.commons.math.distribution.IntegerDistribution;
import org.apache.commons.math.distribution.PascalDistributionImpl;
@ -654,7 +654,7 @@ public class RandomDataImpl implements RandomData, Serializable {
}
/**
* Generates a random value from the {@link FDistributionImpl F Distribution}.
* Generates a random value from the {@link FDistribution F Distribution}.
* This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
* to generate random values.
*
@ -664,12 +664,12 @@ public class RandomDataImpl implements RandomData, Serializable {
* @since 2.2
*/
public double nextF(double numeratorDf, double denominatorDf) {
return nextInversionDeviate(new FDistributionImpl(numeratorDf, denominatorDf));
return nextInversionDeviate(new FDistribution(numeratorDf, denominatorDf));
}
/**
* <p>Generates a random value from the
* {@link org.apache.commons.math.distribution.GammaDistributionImpl Gamma Distribution}.</p>
* {@link org.apache.commons.math.distribution.GammaDistribution Gamma Distribution}.</p>
*
* <p>This implementation uses the following algorithms: </p>
*

3
src/main/java/org/apache/commons/math/stat/inference/OneWayAnovaImpl.java
View File

@ -21,7 +21,6 @@ import java.util.Collection;
import org.apache.commons.math.MathException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.distribution.FDistribution;
import org.apache.commons.math.distribution.FDistributionImpl;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.stat.descriptive.summary.Sum;
import org.apache.commons.math.stat.descriptive.summary.SumOfSquares;
@ -84,7 +83,7 @@ public class OneWayAnovaImpl implements OneWayAnova {
public double anovaPValue(Collection<double[]> categoryData)
throws IllegalArgumentException, MathException {
AnovaStats a = anovaStats(categoryData);
FDistribution fdist = new FDistributionImpl(a.dfbg, a.dfwg);
FDistribution fdist = new FDistribution(a.dfbg, a.dfwg);
return 1.0 - fdist.cumulativeProbability(a.F);
}

16
src/test/java/org/apache/commons/math/distribution/ExponentialDistributionTest.java
View File

@ -43,7 +43,7 @@ public class ExponentialDistributionTest extends ContinuousDistributionAbstractT
/** Creates the default continuous distribution instance to use in tests. */
@Override
public ExponentialDistribution makeDistribution() {
return new ExponentialDistributionImpl(5.0);
return new ExponentialDistribution(5.0);
}
/** Creates the default cumulative probability distribution test input values */
@ -92,14 +92,14 @@ public class ExponentialDistributionTest extends ContinuousDistributionAbstractT
@Test
public void testDensity() {
ExponentialDistribution d1 = new ExponentialDistributionImpl(1);
ExponentialDistribution d1 = new ExponentialDistribution(1);
Assert.assertTrue(Precision.equals(0.0, d1.density(-1e-9), 1));
Assert.assertTrue(Precision.equals(1.0, d1.density(0.0), 1));
Assert.assertTrue(Precision.equals(0.0, d1.density(1000.0), 1));
Assert.assertTrue(Precision.equals(FastMath.exp(-1), d1.density(1.0), 1));
Assert.assertTrue(Precision.equals(FastMath.exp(-2), d1.density(2.0), 1));
ExponentialDistribution d2 = new ExponentialDistributionImpl(3);
ExponentialDistribution d2 = new ExponentialDistribution(3);
Assert.assertTrue(Precision.equals(1/3.0, d2.density(0.0), 1));
// computed using print(dexp(1, rate=1/3), digits=10) in R 2.5
Assert.assertEquals(0.2388437702, d2.density(1.0), 1e-8);
@ -116,19 +116,19 @@ public class ExponentialDistributionTest extends ContinuousDistributionAbstractT
@Test(expected=NotStrictlyPositiveException.class)
public void testPreconditions() {
new ExponentialDistributionImpl(0);
new ExponentialDistribution(0);
}
@Test
public void testMoments() {
final double tol = 1e-9;
ExponentialDistribution dist;
dist = new ExponentialDistributionImpl(11d);
dist = new ExponentialDistribution(11d);
Assert.assertEquals(dist.getNumericalMean(), 11d, tol);
Assert.assertEquals(dist.getNumericalVariance(), 11d * 11d, tol);
dist = new ExponentialDistributionImpl(10.5d);
dist = new ExponentialDistribution(10.5d);
Assert.assertEquals(dist.getNumericalMean(), 10.5d, tol);
Assert.assertEquals(dist.getNumericalVariance(), 10.5d * 10.5d, tol);
}

26
src/test/java/org/apache/commons/math/distribution/FDistributionTest.java
View File

@ -34,7 +34,7 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest {
/** Creates the default continuous distribution instance to use in tests. */
@Override
public FDistribution makeDistribution() {
return new FDistributionImpl(5.0, 6.0);
return new FDistribution(5.0, 6.0);
}
/** Creates the default cumulative probability distribution test input values */
@ -91,13 +91,13 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest {
@Test
public void testPreconditions() {
try {
new FDistributionImpl(0, 1);
new FDistribution(0, 1);
Assert.fail("Expecting NotStrictlyPositiveException for df = 0");
} catch (NotStrictlyPositiveException ex) {
// Expected.
}
try {
new FDistributionImpl(1, 0);
new FDistribution(1, 0);
Assert.fail("Expecting NotStrictlyPositiveException for df = 0");
} catch (NotStrictlyPositiveException ex) {
// Expected.
@ -106,7 +106,7 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest {
@Test
public void testLargeDegreesOfFreedom() throws Exception {
FDistributionImpl fd = new FDistributionImpl(100000, 100000);
FDistribution fd = new FDistribution(100000, 100000);
double p = fd.cumulativeProbability(.999);
double x = fd.inverseCumulativeProbability(p);
Assert.assertEquals(.999, x, 1.0e-5);
@ -114,12 +114,12 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest {
@Test
public void testSmallDegreesOfFreedom() throws Exception {
FDistributionImpl fd = new FDistributionImpl(1, 1);
FDistribution fd = new FDistribution(1, 1);
double p = fd.cumulativeProbability(0.975);
double x = fd.inverseCumulativeProbability(p);
Assert.assertEquals(0.975, x, 1.0e-5);
fd = new FDistributionImpl(1, 2);
fd = new FDistribution(1, 2);
p = fd.cumulativeProbability(0.975);
x = fd.inverseCumulativeProbability(p);
Assert.assertEquals(0.975, x, 1.0e-5);
@ -129,17 +129,17 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest {
public void testMoments() {
final double tol = 1e-9;
FDistribution dist;
dist = new FDistributionImpl(1, 2);
dist = new FDistribution(1, 2);
Assert.assertTrue(Double.isNaN(dist.getNumericalMean()));
Assert.assertTrue(Double.isNaN(dist.getNumericalVariance()));
dist = new FDistributionImpl(1, 3);
dist = new FDistribution(1, 3);
Assert.assertEquals(dist.getNumericalMean(), 3d / (3d - 2d), tol);
Assert.assertTrue(Double.isNaN(dist.getNumericalVariance()));
dist = new FDistributionImpl(1, 5);
dist = new FDistribution(1, 5);
Assert.assertEquals(dist.getNumericalMean(), 5d / (5d - 2d), tol);
Assert.assertEquals(dist.getNumericalVariance(), (2d * 5d * 5d * 4d) / 9d, tol);
Assert.assertEquals(dist.getNumericalVariance(), (2d * 5d * 5d * 4d) / 9d, tol);
}
}

22
src/test/java/org/apache/commons/math/distribution/GammaDistributionTest.java
View File

@ -35,7 +35,7 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest {
/** Creates the default continuous distribution instance to use in tests. */
@Override
public GammaDistribution makeDistribution() {
return new GammaDistributionImpl(4d, 2d);
return new GammaDistribution(4d, 2d);
}
/** Creates the default cumulative probability distribution test input values */
@ -77,13 +77,13 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest {
@Test
public void testPreconditions() {
try {
new GammaDistributionImpl(0, 1);
new GammaDistribution(0, 1);
Assert.fail("Expecting NotStrictlyPositiveException for alpha = 0");
} catch (NotStrictlyPositiveException ex) {
// Expected.
}
try {
new GammaDistributionImpl(1, 0);
new GammaDistribution(1, 0);
Assert.fail("Expecting NotStrictlyPositiveException for alpha = 0");
} catch (NotStrictlyPositiveException ex) {
// Expected.
@ -108,13 +108,13 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest {
}
private void testProbability(double x, double a, double b, double expected) throws Exception {
GammaDistribution distribution = new GammaDistributionImpl( a, b );
GammaDistribution distribution = new GammaDistribution( a, b );
double actual = distribution.cumulativeProbability(x);
Assert.assertEquals("probability for " + x, expected, actual, 10e-4);
}
private void testValue(double expected, double a, double b, double p) throws Exception {
GammaDistribution distribution = new GammaDistributionImpl( a, b );
GammaDistribution distribution = new GammaDistribution( a, b );
double actual = distribution.inverseCumulativeProbability(p);
Assert.assertEquals("critical value for " + p, expected, actual, 10e-4);
}
@ -141,7 +141,7 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest {
}
private void checkDensity(double alpha, double rate, double[] x, double[] expected) {
GammaDistribution d = new GammaDistributionImpl(alpha, 1 / rate);
GammaDistribution d = new GammaDistribution(alpha, 1 / rate);
for (int i = 0; i < x.length; i++) {
Assert.assertEquals(expected[i], d.density(x[i]), 1e-5);
}
@ -158,12 +158,12 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest {
public void testMoments() {
final double tol = 1e-9;
GammaDistribution dist;
dist = new GammaDistributionImpl(1, 2);
dist = new GammaDistribution(1, 2);
Assert.assertEquals(dist.getNumericalMean(), 2, tol);
Assert.assertEquals(dist.getNumericalVariance(), 4, tol);
dist = new GammaDistributionImpl(1.1, 4.2);
Assert.assertEquals(dist.getNumericalVariance(), 4, tol);
dist = new GammaDistribution(1.1, 4.2);
Assert.assertEquals(dist.getNumericalMean(), 1.1d * 4.2d, tol);
Assert.assertEquals(dist.getNumericalVariance(), 1.1d * 4.2d * 4.2d, tol);
}

16
src/test/java/org/apache/commons/math/random/RandomDataTest.java
View File

@ -30,9 +30,9 @@ import org.apache.commons.math.distribution.BinomialDistribution;
import org.apache.commons.math.distribution.BinomialDistributionTest;
import org.apache.commons.math.distribution.CauchyDistribution;
import org.apache.commons.math.distribution.ChiSquaredDistribution;
import org.apache.commons.math.distribution.ExponentialDistributionImpl;
import org.apache.commons.math.distribution.FDistributionImpl;
import org.apache.commons.math.distribution.GammaDistributionImpl;
import org.apache.commons.math.distribution.ExponentialDistribution;
import org.apache.commons.math.distribution.FDistribution;
import org.apache.commons.math.distribution.GammaDistribution;
import org.apache.commons.math.distribution.HypergeometricDistributionImpl;
import org.apache.commons.math.distribution.HypergeometricDistributionTest;
import org.apache.commons.math.distribution.PascalDistributionImpl;
@ -616,7 +616,7 @@ public class RandomDataTest {
long[] counts;
// Mean 1
quartiles = TestUtils.getDistributionQuartiles(new ExponentialDistributionImpl(1));
quartiles = TestUtils.getDistributionQuartiles(new ExponentialDistribution(1));
counts = new long[4];
randomData.reSeed(1000);
for (int i = 0; i < 1000; i++) {
@ -626,7 +626,7 @@ public class RandomDataTest {
TestUtils.assertChiSquareAccept(expected, counts, 0.001);
// Mean 5
quartiles = TestUtils.getDistributionQuartiles(new ExponentialDistributionImpl(5));
quartiles = TestUtils.getDistributionQuartiles(new ExponentialDistribution(5));
counts = new long[4];
randomData.reSeed(1000);
for (int i = 0; i < 1000; i++) {
@ -896,7 +896,7 @@ public class RandomDataTest {
@Test
public void testNextF() throws Exception {
double[] quartiles = TestUtils.getDistributionQuartiles(new FDistributionImpl(12, 5));
double[] quartiles = TestUtils.getDistributionQuartiles(new FDistribution(12, 5));
long[] counts = new long[4];
randomData.reSeed(1000);
for (int i = 0; i < 1000; i++) {
@ -912,7 +912,7 @@ public class RandomDataTest {
long[] counts;
// Tests shape > 1, one case in the rejection sampling
quartiles = TestUtils.getDistributionQuartiles(new GammaDistributionImpl(4, 2));
quartiles = TestUtils.getDistributionQuartiles(new GammaDistribution(4, 2));
counts = new long[4];
randomData.reSeed(1000);
for (int i = 0; i < 1000; i++) {
@ -922,7 +922,7 @@ public class RandomDataTest {
TestUtils.assertChiSquareAccept(expected, counts, 0.001);
// Tests shape <= 1, another case in the rejection sampling
quartiles = TestUtils.getDistributionQuartiles(new GammaDistributionImpl(0.3, 3));
quartiles = TestUtils.getDistributionQuartiles(new GammaDistribution(0.3, 3));
counts = new long[4];
randomData.reSeed(1000);
for (int i = 0; i < 1000; i++) {