Hi
I've been looking around for open source mathematical statistics software in
Java in the last couple of weeks and have tested the commons-math package.
After looking into it I have a few questions and comments.
1. The MathUtils.factorial(final int n) method throws an
IllegalArgumentException for n=0 which is clearly incorrect since 0!=1 by
definition.
2. I think the distinction between discrete and continuous probability
distributions manifested in the interfaces DiscreteDistribution and
ContinuousDistribution is rather questionable There is no requirement that
a discrete distribution only takes on integer values so the methods of
interface DiscreteDistribution doesn't cover all discrete distributions. On
the other hand, all of the methods of ContinuousDistribution holds equally
well for a discrete probability distribution. In my opinion, a more
appropriate approach would be to rename ContinuousDistribution to
ProbabilityDistribution and drop the DiscreteDistribution interface. Of
course, it could be practical to have convenience methods that takes integer
arguments for the probability densities for certain distributions but then
you can define a new interface IntegerValuedDistribution like
public interface IntegerValuedDistribution extends ProbabilityDistribution {
double probability(int x);
double cumulativeProbability(int x) throws MathException;
}
3. Wouldn't it be nice to have convenience methods for calculating
the moments for each distribution? Something along the lines of
public double getMoment(int order) throws MomentDoesNotExistException;
for calculating the moments EX^order (if they exist)?
At least there could be methods for obtaining the mean and variance of a
distribution.
4. Since the chi-squared and exponential distributions are just special cases
of the gamma distribution, there is no need to have separate implementation
classes for these. In my opinion, one should avoid having multiple
implementations of the same distribution unless there is some strong reason
for it.
5. There are quite a lot of elementary distributions missing.
I wrote an implementation of the Poisson distribution while testing the
package and have attached the files for it.
Best regards,
Frank N
/*
* Copyright 2003-2004 The Apache Software Foundation.
*
* Licensed 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 org.apache.commons.math.MathException;
/**
* <code>The Poisson distribution</code>
*
* <p>
* References:
* <ul>
* <li><a href="http://mathworld.wolfram.com/PoissonDistribution.html";>
* Poisson distribution</a></li>
* </ul>
* </p>
*
* @version $Id: $
*/
public interface PoissonDistribution extends DiscreteDistribution {
/**
* Get the Poisson parameter for the distribution.
* @return the Poisson parameter for the distribution.
*/
public double getParameter();
/**
* Set the Poisson parameter for the distribution.
*/
public void setParameter(double p);
/**
* Calculates the Poisson distribution function using a normal approximation.
* @param x the upper bound, inclusive
* @return the distribution function value calculated using a normal approximation
* @throws MathException
*/
public double normalApproximateProbability(int x) throws MathException;
}/*
* Copyright 2003-2004 The Apache Software Foundation.
*
* Licensed 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 java.text.DecimalFormat;
import org.apache.commons.math.MathException;
import org.apache.commons.math.util.MathUtils;
/**
* <code>PoissonDistributionImpl</code>
*
* @version $Id: $
*/
public class PoissonDistributionImpl extends AbstractDiscreteDistribution
implements PoissonDistribution, Serializable {
/**
* Holds the Poisson parameter for the distribution.
*/
private double parameter;
/**
* Create a new Poisson distribution with the given the parameter.
* @param the Poisson parameter
*/
public PoissonDistributionImpl(double p) {
super();
setParameter(p);
}
/**
* Get the Poisson parameter for the distribution.
* @return the Poisson parameter for the distribution.
*/
public double getParameter() {
return this.parameter;
}
/**
* Set the Poisson parameter for the distribution.
*/
public void setParameter(double p) {
if (p <= 0) {
throw new IllegalArgumentException(
"The Poisson parameter must be positive");
}
this.parameter = p;
}
/**
* The probability density function P(X = x) for a Poisson distribution.
*/
public double probability(int x) {
if (x < 0 || x == Integer.MAX_VALUE) {
return 0;
}
return Math.pow(getParameter(), x)
/ (x > 0 ? MathUtils.factorial(x) : 1) * Math.exp(-parameter);
}
/**
* The probability distribution function P(X <= x) for a Poisson distribution.
*/
public double cumulativeProbability(int x) throws MathException {
if (x < 0) {
return 0;
}
if (x == Integer.MAX_VALUE) {
return 1;
}
double result = 0;
for (int i = 0; i <= x; i++) {
result += probability(i);
}
return result;
}
/**
* For this distribution, X, this method returns the largest x, such
* that P(X ≤ x) ≤ <code>p</code>.
*
* @param p the desired probability
* @return the largest x such that P(X ≤ x) <= p
* @throws MathException if the inverse cumulative probability can not be
* computed due to convergence or other numerical errors.
* @throws IllegalArgumentException if p < 0 or p > 1
*/
public int inverseCumulativeProbability(final double p)
throws MathException {
if (p < 0.0 || p > 1.0) {
throw new IllegalArgumentException(
"p must be between 0 and 1.0 (inclusive)");
}
if (p == 1d) {
throw new MathException(
"No upper bound exists for an inverse cumulative probability of 1.0 for a Poisson distribution.");
}
// check for degenerate case p=0
if (p == 0) {
return -1;
}
// use a DecimalFormat to represent the precision that should be used
// in the calculations
DecimalFormat numFormat = createFormat(p);
int x = -1;
double pm = 0;
do {
x++;
pm = cumulativeProbability(x);
// convert to the number format of p
pm = new Double(numFormat.format(pm)).doubleValue();
} while (pm < p);
int x0;
if (pm <= p) {
x0 = x;
}
else {
x0 = x - 1;
}
return x0;
}
/**
* Calculates the Poisson distribution function using a normal approximation.
* @param x the upper bound, inclusive
* @return the distribution function value calculated using a normal approximation
* @throws MathException
*/
public double normalApproximateProbability(int x) throws MathException {
NormalDistribution normal = DistributionFactory.newInstance()
.createNormalDistribution(getParameter(),
Math.sqrt(getParameter()));
// calculate the probability using half-correction
return normal.cumulativeProbability(x + 0.5);
}
/**
* Access the domain value lower bound, based on <code>p</code>, used to
* bracket a CDF root. This method is used by
* [EMAIL PROTECTED] #inverseCumulativeProbability(double)} to find critical values.
*
* @param p the desired probability for the critical value
*/
protected int getDomainLowerBound(double p) {
return 0;
}
/**
* Access the domain value upper bound, based on <code>p</code>, used to
* bracket a CDF root. This method is used by
* [EMAIL PROTECTED] #inverseCumulativeProbability(double)} to find critical values.
*
* @param p the desired probability for the critical value
*/
protected int getDomainUpperBound(double p) {
return Integer.MAX_VALUE;
}
/**
* This method creates a <code>DecimalFormat</code> corresponding to the
* decimal precision of a given double.
* @param d
*/
protected DecimalFormat createFormat(double d) {
String s = new Double(d).toString().replaceAll("[1-9]", "0");
return new DecimalFormat(s);
}
}/*
* Copyright 2003-2004 The Apache Software Foundation.
*
* Licensed 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.text.DecimalFormat;
import org.apache.commons.math.MathException;
/**
* <code>PoissonDistributionTest</code>
*
* @version $Id: $
*/
public class PoissonDistributionTest extends DiscreteDistributionAbstractTest {
/**
* Poisson parameter value for the test distribution.
*/
private static final double DEFAULT_TEST_POISSON_PARAMETER = 4.0;
/**
* Constructor.
* @param name
*/
public PoissonDistributionTest(String name) {
super(name);
}
/**
* Creates the default discrete distribution instance to use in tests.
*/
public DiscreteDistribution makeDistribution() {
return new PoissonDistributionImpl(DEFAULT_TEST_POISSON_PARAMETER);
}
/**
* Creates the default probability density test input values.
*/
public int[] makeDensityTestPoints() {
return new int[] { -1, 0, 1, 2, 3, 4, 5, 10 };
}
/**
* Creates the default probability density test expected values.
*/
public double[] makeDensityTestValues() {
return new double[] { 0d, 0.01832d, 0.07326d, 0.14653d, 0.19534d,
0.19534d, 0.15629d, 0.0053d };
}
/**
* Creates the default cumulative probability density test input values.
*/
public int[] makeCumulativeTestPoints() {
return new int[] { -1, 0, 1, 2, 3, 4, 5, 10 };
}
/**
* Creates the default cumulative probability density test expected values.
*/
public double[] makeCumulativeTestValues() {
return new double[] { 0d, 0.01832d, 0.09158d, 0.23810d, 0.43347d,
0.62884d, 0.78513d, 0.9972d };
}
/**
* Creates the default inverse cumulative probability test input values.
*/
public double[] makeInverseCumulativeTestPoints() {
return new double[] { 0.0d, 0.01832d, 0.09158d, 0.23810d, 0.43347d,
0.62884d, 0.78513d };
}
/**
* Creates the default inverse cumulative probability density test expected values.
*/
public int[] makeInverseCumulativeTestValues() {
return new int[] { -1, 0, 1, 2, 3, 4, 5 };
}
/**
* Test the normal approximation of the Poisson distribution by
* calculating P(90 ≤ X ≤ 110) if X is Po(100).
*/
public void testNormalApproximateProbability() throws Exception {
PoissonDistribution dist = new PoissonDistributionImpl(100);
double result = dist.normalApproximateProbability(110)
- dist.normalApproximateProbability(89);
String strResult = new DecimalFormat("0.000").format(result);
assertEquals("0.706", strResult);
}
/**
* Test the degenerate case of a 1.0 inverse cumulative probability.
* @throws Exception
*/
public void testDegenerateInverseCumulativeProbability() throws Exception {
try {
PoissonDistribution dist = new PoissonDistributionImpl(
DEFAULT_TEST_POISSON_PARAMETER);
dist.inverseCumulativeProbability(1.0d);
fail("Degenerate Poisson inverse cumulative probability test should have failed.");
}
catch (MathException me) {
assertTrue(true);
}
}
}
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
