commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From brentwor...@apache.org
Subject svn commit: r830745 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/distribution/ site/xdoc/ test/java/org/apache/commons/math/distribution/
Date Wed, 28 Oct 2009 19:59:22 GMT
Author: brentworden
Date: Wed Oct 28 19:59:21 2009
New Revision: 830745

URL: http://svn.apache.org/viewvc?rev=830745&view=rev
Log:
MATH-311.  Changed probability calculations for Binomial, Poisson, and Hypergeometric distributions
to use Catherine Loader's saddle point approximations

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/SaddlePointExpansion.java
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/BinomialDistributionImpl.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/distribution/BinomialDistributionTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/BinomialDistributionImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/BinomialDistributionImpl.java?rev=830745&r1=830744&r2=830745&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/BinomialDistributionImpl.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/BinomialDistributionImpl.java
Wed Oct 28 19:59:21 2009
@@ -1,18 +1,15 @@
 /*
  * 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.
+ * 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;
 
@@ -25,12 +22,12 @@
 
 /**
  * The default implementation of {@link BinomialDistribution}.
- *
- * @version $Revision$ $Date$
+ * 
+ * @version $Revision$ $Date: 2009-09-05 12:36:48 -0500 (Sat, 05 Sep
+ *          2009) $
  */
-public class BinomialDistributionImpl
-    extends AbstractIntegerDistribution
-    implements BinomialDistribution, Serializable {
+public class BinomialDistributionImpl extends AbstractIntegerDistribution
+        implements BinomialDistribution, Serializable {
 
     /** Serializable version identifier */
     private static final long serialVersionUID = 6751309484392813623L;
@@ -44,6 +41,7 @@
     /**
      * Create a binomial distribution with the given number of trials and
      * probability of success.
+     * 
      * @param trials the number of trials.
      * @param p the probability of success.
      */
@@ -55,6 +53,7 @@
 
     /**
      * Access the number of trials for this distribution.
+     * 
      * @return the number of trials.
      */
     public int getNumberOfTrials() {
@@ -63,6 +62,7 @@
 
     /**
      * Access the probability of success for this distribution.
+     * 
      * @return the probability of success.
      */
     public double getProbabilityOfSuccess() {
@@ -71,28 +71,30 @@
 
     /**
      * Change the number of trials for this distribution.
+     * 
      * @param trials the new number of trials.
      * @throws IllegalArgumentException if <code>trials</code> is not a valid
-     *         number of trials.
+     *             number of trials.
      */
     public void setNumberOfTrials(int trials) {
         if (trials < 0) {
             throw MathRuntimeException.createIllegalArgumentException(
-                  "number of trials must be non-negative ({0})", trials);
+                    "number of trials must be non-negative ({0})", trials);
         }
         numberOfTrials = trials;
     }
 
     /**
      * Change the probability of success for this distribution.
+     * 
      * @param p the new probability of success.
      * @throws IllegalArgumentException if <code>p</code> is not a valid
-     *         probability.
+     *             probability.
      */
     public void setProbabilityOfSuccess(double p) {
         if (p < 0.0 || p > 1.0) {
             throw MathRuntimeException.createIllegalArgumentException(
-                  "{0} out of [{1}, {2}] range", p, 0.0, 1.0);
+                    "{0} out of [{1}, {2}] range", p, 0.0, 1.0);
         }
         probabilityOfSuccess = p;
     }
@@ -100,10 +102,10 @@
     /**
      * Access the domain value lower bound, based on <code>p</code>, used to
      * bracket a PDF root.
-     *
+     * 
      * @param p the desired probability for the critical value
-     * @return domain value lower bound, i.e.
-     *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
+     * @return domain value lower bound, i.e. P(X &lt; <i>lower bound</i>)
&lt;
+     *         <code>p</code>
      */
     @Override
     protected int getDomainLowerBound(double p) {
@@ -113,10 +115,10 @@
     /**
      * Access the domain value upper bound, based on <code>p</code>, used to
      * bracket a PDF root.
-     *
+     * 
      * @param p the desired probability for the critical value
-     * @return domain value upper bound, i.e.
-     *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
+     * @return domain value upper bound, i.e. P(X &lt; <i>upper bound</i>)
&gt;
+     *         <code>p</code>
      */
     @Override
     protected int getDomainUpperBound(double p) {
@@ -125,10 +127,11 @@
 
     /**
      * For this distribution, X, this method returns P(X &le; x).
+     * 
      * @param x the value at which the PDF is evaluated.
      * @return PDF for this distribution.
-     * @throws MathException if the cumulative probability can not be
-     *            computed due to convergence or other numerical errors.
+     * @throws MathException if the cumulative probability can not be computed
+     *             due to convergence or other numerical errors.
      */
     @Override
     public double cumulativeProbability(int x) throws MathException {
@@ -138,18 +141,15 @@
         } else if (x >= getNumberOfTrials()) {
             ret = 1.0;
         } else {
-            ret =
-                1.0 - Beta.regularizedBeta(
-                        getProbabilityOfSuccess(),
-                        x + 1.0,
-                        getNumberOfTrials() - x);
+            ret = 1.0 - Beta.regularizedBeta(getProbabilityOfSuccess(),
+                    x + 1.0, getNumberOfTrials() - x);
         }
         return ret;
     }
 
     /**
      * For this distribution, X, this method returns P(X = x).
-     *
+     * 
      * @param x the value at which the PMF is evaluated.
      * @return PMF for this distribution.
      */
@@ -158,30 +158,30 @@
         if (x < 0 || x > getNumberOfTrials()) {
             ret = 0.0;
         } else {
-            ret = MathUtils.binomialCoefficientDouble(
-                    getNumberOfTrials(), x) *
-                  Math.pow(getProbabilityOfSuccess(), x) *
-                  Math.pow(1.0 - getProbabilityOfSuccess(),
-                        getNumberOfTrials() - x);
+            ret = Math.exp(SaddlePointExpansion.logBinomialProbability(x,
+                    numberOfTrials, probabilityOfSuccess,
+                    1.0 - probabilityOfSuccess));
         }
         return ret;
     }
 
     /**
-     * For this distribution, X, this method returns the largest x, such
-     * that P(X &le; x) &le; <code>p</code>.
+     * For this distribution, X, this method returns the largest x, such that
+     * P(X &le; x) &le; <code>p</code>.
      * <p>
      * Returns <code>-1</code> for p=0 and <code>Integer.MAX_VALUE</code>
for
-     * p=1.</p>
-     *
+     * p=1.
+     * </p>
+     * 
      * @param p the desired probability
      * @return the largest x such that P(X &le; x) <= p
      * @throws MathException if the inverse cumulative probability can not be
-     *            computed due to convergence or other numerical errors.
+     *             computed due to convergence or other numerical errors.
      * @throws IllegalArgumentException if p < 0 or p > 1
      */
     @Override
-    public int inverseCumulativeProbability(final double p) throws MathException {
+    public int inverseCumulativeProbability(final double p)
+            throws MathException {
         // handle extreme values explicitly
         if (p == 0) {
             return -1;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java?rev=830745&r1=830744&r2=830745&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java
Wed Oct 28 19:59:21 2009
@@ -1,18 +1,15 @@
 /*
  * 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.
+ * 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;
@@ -24,12 +21,12 @@
 
 /**
  * The default implementation of {@link HypergeometricDistribution}.
- *
- * @version $Revision$ $Date$
+ * 
+ * @version $Revision$ $Date: 2009-09-05 12:36:48 -0500 (Sat, 05 Sep
+ *          2009) $
  */
 public class HypergeometricDistributionImpl extends AbstractIntegerDistribution
-    implements HypergeometricDistribution, Serializable
-{
+        implements HypergeometricDistribution, Serializable {
 
     /** Serializable version identifier */
     private static final long serialVersionUID = -436928820673516179L;
@@ -46,22 +43,25 @@
     /**
      * Construct a new hypergeometric distribution with the given the population
      * size, the number of successes in the population, and the sample size.
+     * 
      * @param populationSize the population size.
      * @param numberOfSuccesses number of successes in the population.
      * @param sampleSize the sample size.
      */
     public HypergeometricDistributionImpl(int populationSize,
-        int numberOfSuccesses, int sampleSize) {
+            int numberOfSuccesses, int sampleSize) {
         super();
         if (numberOfSuccesses > populationSize) {
-            throw MathRuntimeException.createIllegalArgumentException(
-                "number of successes ({0}) must be less than or equal to population size
({1})",
-                numberOfSuccesses, populationSize);
+            throw MathRuntimeException
+                    .createIllegalArgumentException(
+                            "number of successes ({0}) must be less than or equal to population
size ({1})",
+                            numberOfSuccesses, populationSize);
         }
         if (sampleSize > populationSize) {
-            throw MathRuntimeException.createIllegalArgumentException(
-                  "sample size ({0}) must be less than or equal to population size ({1})",
-                  sampleSize, populationSize);
+            throw MathRuntimeException
+                    .createIllegalArgumentException(
+                            "sample size ({0}) must be less than or equal to population size
({1})",
+                            sampleSize, populationSize);
         }
         setPopulationSize(populationSize);
         setSampleSize(sampleSize);
@@ -70,6 +70,7 @@
 
     /**
      * For this distribution, X, this method returns P(X &le; x).
+     * 
      * @param x the value at which the PDF is evaluated.
      * @return PDF for this distribution.
      */
@@ -84,7 +85,7 @@
         int[] domain = getDomain(n, m, k);
         if (x < domain[0]) {
             ret = 0.0;
-        } else if(x >= domain[1]) {
+        } else if (x >= domain[1]) {
             ret = 1.0;
         } else {
             ret = innerCumulativeProbability(domain[0], x, 1, n, m, k);
@@ -95,40 +96,38 @@
 
     /**
      * Return the domain for the given hypergeometric distribution parameters.
+     * 
      * @param n the population size.
      * @param m number of successes in the population.
      * @param k the sample size.
      * @return a two element array containing the lower and upper bounds of the
      *         hypergeometric distribution.
      */
-    private int[] getDomain(int n, int m, int k){
-        return new int[]{
-            getLowerDomain(n, m, k),
-            getUpperDomain(m, k)
-        };
+    private int[] getDomain(int n, int m, int k) {
+        return new int[] { getLowerDomain(n, m, k), getUpperDomain(m, k) };
     }
 
     /**
      * Access the domain value lower bound, based on <code>p</code>, used to
      * bracket a PDF root.
-     *
+     * 
      * @param p the desired probability for the critical value
-     * @return domain value lower bound, i.e.
-     *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
+     * @return domain value lower bound, i.e. P(X &lt; <i>lower bound</i>)
&lt;
+     *         <code>p</code>
      */
     @Override
     protected int getDomainLowerBound(double p) {
         return getLowerDomain(getPopulationSize(), getNumberOfSuccesses(),
-            getSampleSize());
+                getSampleSize());
     }
 
     /**
      * Access the domain value upper bound, based on <code>p</code>, used to
      * bracket a PDF root.
-     *
+     * 
      * @param p the desired probability for the critical value
-     * @return domain value upper bound, i.e.
-     *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
+     * @return domain value upper bound, i.e. P(X &lt; <i>upper bound</i>)
&gt;
+     *         <code>p</code>
      */
     @Override
     protected int getDomainUpperBound(double p) {
@@ -138,6 +137,7 @@
     /**
      * Return the lowest domain value for the given hypergeometric distribution
      * parameters.
+     * 
      * @param n the population size.
      * @param m number of successes in the population.
      * @param k the sample size.
@@ -149,6 +149,7 @@
 
     /**
      * Access the number of successes.
+     * 
      * @return the number of successes.
      */
     public int getNumberOfSuccesses() {
@@ -157,6 +158,7 @@
 
     /**
      * Access the population size.
+     * 
      * @return the population size.
      */
     public int getPopulationSize() {
@@ -165,6 +167,7 @@
 
     /**
      * Access the sample size.
+     * 
      * @return the sample size.
      */
     public int getSampleSize() {
@@ -174,32 +177,42 @@
     /**
      * Return the highest domain value for the given hypergeometric distribution
      * parameters.
+     * 
      * @param m number of successes in the population.
      * @param k the sample size.
      * @return the highest domain value of the hypergeometric distribution.
      */
-    private int getUpperDomain(int m, int k){
+    private int getUpperDomain(int m, int k) {
         return Math.min(k, m);
     }
 
     /**
      * For this distribution, X, this method returns P(X = x).
-     *
+     * 
      * @param x the value at which the PMF is evaluated.
      * @return PMF for this distribution.
      */
     public double probability(int x) {
         double ret;
 
-        int n = getPopulationSize();
-        int m = getNumberOfSuccesses();
+        int m = getPopulationSize();
+        int s = getNumberOfSuccesses();
+        int f = m - s;
         int k = getSampleSize();
 
-        int[] domain = getDomain(n, m, k);
-        if(x < domain[0] || x > domain[1]){
+        int[] domain = getDomain(m, s, k);
+        if (x < domain[0] || x > domain[1]) {
             ret = 0.0;
         } else {
-            ret = probability(n, m, k, x);
+            double p = (double) sampleSize / (double) m;
+            double q = (double) (m - sampleSize) / (double) m;
+            double p1 = SaddlePointExpansion.logBinomialProbability(x,
+                    numberOfSuccesses, p, q);
+            double p2 = SaddlePointExpansion.logBinomialProbability(sampleSize
+                    - x, f, p, q);
+            double p3 = SaddlePointExpansion.logBinomialProbability(sampleSize,
+                    m, p, q);
+            ret = Math.exp(p1 + p2 - p3);
         }
 
         return ret;
@@ -208,7 +221,7 @@
     /**
      * For the distribution, X, defined by the given hypergeometric distribution
      * parameters, this method returns P(X = x).
-     *
+     * 
      * @param n the population size.
      * @param m number of successes in the population.
      * @param k the sample size.
@@ -216,55 +229,56 @@
      * @return PMF for the distribution.
      */
     private double probability(int n, int m, int k, int x) {
-        return Math.exp(MathUtils.binomialCoefficientLog(m, x) +
-            MathUtils.binomialCoefficientLog(n - m, k - x) -
-            MathUtils.binomialCoefficientLog(n, k));
+        return Math.exp(MathUtils.binomialCoefficientLog(m, x)
+                + MathUtils.binomialCoefficientLog(n - m, k - x)
+                - MathUtils.binomialCoefficientLog(n, k));
     }
 
     /**
      * Modify the number of successes.
+     * 
      * @param num the new number of successes.
      * @throws IllegalArgumentException if <code>num</code> is negative.
      */
     public void setNumberOfSuccesses(int num) {
-        if(num < 0){
+        if (num < 0) {
             throw MathRuntimeException.createIllegalArgumentException(
-                  "number of successes must be non-negative ({0})",
-                  num);
+                    "number of successes must be non-negative ({0})", num);
         }
         numberOfSuccesses = num;
     }
 
     /**
      * Modify the population size.
+     * 
      * @param size the new population size.
      * @throws IllegalArgumentException if <code>size</code> is not positive.
      */
     public void setPopulationSize(int size) {
-        if(size <= 0){
+        if (size <= 0) {
             throw MathRuntimeException.createIllegalArgumentException(
-                  "population size must be positive ({0})",
-                  size);
+                    "population size must be positive ({0})", size);
         }
         populationSize = size;
     }
 
     /**
      * Modify the sample size.
+     * 
      * @param size the new sample size.
      * @throws IllegalArgumentException if <code>size</code> is negative.
      */
     public void setSampleSize(int size) {
         if (size < 0) {
             throw MathRuntimeException.createIllegalArgumentException(
-                  "sample size must be positive ({0})",
-                  size);
+                    "sample size must be positive ({0})", size);
         }
         sampleSize = size;
     }
 
     /**
      * For this distribution, X, this method returns P(X &ge; x).
+     * 
      * @param x the value at which the CDF is evaluated.
      * @return upper tail CDF for this distribution.
      * @since 1.1
@@ -279,7 +293,7 @@
         int[] domain = getDomain(n, m, k);
         if (x < domain[0]) {
             ret = 1.0;
-        } else if(x > domain[1]) {
+        } else if (x > domain[1]) {
             ret = 0.0;
         } else {
             ret = innerCumulativeProbability(domain[1], x, -1, n, m, k);
@@ -289,21 +303,21 @@
     }
 
     /**
-     * For this distribution, X, this method returns P(x0 &le; X &le; x1).  This
+     * For this distribution, X, this method returns P(x0 &le; X &le; x1). This
      * probability is computed by summing the point probabilities for the values
      * x0, x0 + 1, x0 + 2, ..., x1, in the order directed by dx.
+     * 
      * @param x0 the inclusive, lower bound
      * @param x1 the inclusive, upper bound
      * @param dx the direction of summation. 1 indicates summing from x0 to x1.
-     *           0 indicates summing from x1 to x0.
+     *            0 indicates summing from x1 to x0.
      * @param n the population size.
      * @param m number of successes in the population.
      * @param k the sample size.
      * @return P(x0 &le; X &le; x1).
      */
-    private double innerCumulativeProbability(
-        int x0, int x1, int dx, int n, int m, int k)
-    {
+    private double innerCumulativeProbability(int x0, int x1, int dx, int n,
+            int m, int k) {
         double ret = probability(n, m, k, x0);
         while (x0 != x1) {
             x0 += dx;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java?rev=830745&r1=830744&r2=830745&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java
Wed Oct 28 19:59:21 2009
@@ -1,18 +1,15 @@
 /*
  * 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.
+ * 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;
 
@@ -25,8 +22,9 @@
 
 /**
  * Implementation for the {@link PoissonDistribution}.
- *
- * @version $Revision$ $Date$
+ * 
+ * @version $Revision$ $Date: 2009-09-05 12:36:48 -0500 (Sat, 05 Sep
+ *          2009) $
  */
 public class PoissonDistributionImpl extends AbstractIntegerDistribution
         implements PoissonDistribution, Serializable {
@@ -43,10 +41,9 @@
     private double mean;
 
     /**
-     * Create a new Poisson distribution with the given the mean.
-     * The mean value must be positive; otherwise an
-     * <code>IllegalArgument</code> is thrown.
-     *
+     * Create a new Poisson distribution with the given the mean. The mean value
+     * must be positive; otherwise an <code>IllegalArgument</code> is thrown.
+     * 
      * @param p the Poisson mean
      * @throws IllegalArgumentException if p &le; 0
      */
@@ -55,10 +52,9 @@
     }
 
     /**
-     * Create a new Poisson distribution with the given the mean.
-     * The mean value must be positive; otherwise an
-     * <code>IllegalArgument</code> is thrown.
-     *
+     * Create a new Poisson distribution with the given the mean. The mean value
+     * must be positive; otherwise an <code>IllegalArgument</code> is thrown.
+     * 
      * @param p the Poisson mean
      * @param z a normal distribution used to compute normal approximations.
      * @throws IllegalArgumentException if p &le; 0
@@ -72,7 +68,7 @@
 
     /**
      * Get the Poisson mean for the distribution.
-     *
+     * 
      * @return the Poisson mean for the distribution.
      */
     public double getMean() {
@@ -80,18 +76,16 @@
     }
 
     /**
-     * Set the Poisson mean for the distribution.
-     * The mean value must be positive; otherwise an
-     * <code>IllegalArgument</code> is thrown.
-     *
+     * Set the Poisson mean for the distribution. The mean value must be
+     * positive; otherwise an <code>IllegalArgument</code> is thrown.
+     * 
      * @param p the Poisson mean value
      * @throws IllegalArgumentException if p &le; 0
      */
     public void setMean(double p) {
         if (p <= 0) {
             throw MathRuntimeException.createIllegalArgumentException(
-                  "the Poisson mean must be positive ({0})",
-                  p);
+                    "the Poisson mean must be positive ({0})", p);
         }
         this.mean = p;
         normal.setMean(p);
@@ -100,25 +94,34 @@
 
     /**
      * The probability mass function P(X = x) for a Poisson distribution.
-     *
-     * @param x the value at which the probability density function is evaluated.
+     * 
+     * @param x the value at which the probability density function is
+     *            evaluated.
      * @return the value of the probability mass function at x
      */
     public double probability(int x) {
+        double ret;
         if (x < 0 || x == Integer.MAX_VALUE) {
-            return 0;
+            ret = 0.0;
+        } else if (x == 0) {
+            ret = Math.exp(-mean);
+        } else {
+            ret = Math.exp(-SaddlePointExpansion.getStirlingError(x)
+                    - SaddlePointExpansion.getDeviancePart(x, mean))
+                    / Math.sqrt(2.0 * Math.PI * x); // TODO make MathUtils.PI
+                                                    // public
         }
-        return Math.pow(getMean(), x) /
-            MathUtils.factorialDouble(x) * Math.exp(-mean);
+        return ret;
     }
 
     /**
-     * The probability distribution function P(X <= x) for a Poisson distribution.
-     *
+     * The probability distribution function P(X <= x) for a Poisson
+     * distribution.
+     * 
      * @param x the value at which the PDF is evaluated.
      * @return Poisson distribution function evaluated at x
-     * @throws MathException if the cumulative probability can not be
-     *            computed due to convergence or other numerical errors.
+     * @throws MathException if the cumulative probability can not be computed
+     *             due to convergence or other numerical errors.
      */
     @Override
     public double cumulativeProbability(int x) throws MathException {
@@ -128,21 +131,24 @@
         if (x == Integer.MAX_VALUE) {
             return 1;
         }
-        return Gamma.regularizedGammaQ((double)x + 1, mean,
-                1E-12, Integer.MAX_VALUE);
+        return Gamma.regularizedGammaQ((double) x + 1, mean, 1E-12,
+                Integer.MAX_VALUE);
     }
 
     /**
      * Calculates the Poisson distribution function using a normal
-     * approximation.  The <code>N(mean, sqrt(mean))</code>
-     * distribution is used to approximate the Poisson distribution.
+     * approximation. The <code>N(mean, sqrt(mean))</code> distribution is used
+     * to approximate the Poisson distribution.
      * <p>
      * The computation uses "half-correction" -- evaluating the normal
-     * distribution function at <code>x + 0.5</code></p>
-     *
+     * distribution function at <code>x + 0.5</code>
+     * </p>
+     * 
      * @param x the upper bound, inclusive
-     * @return the distribution function value calculated using a normal approximation
-     * @throws MathException if an error occurs computing the normal approximation
+     * @return the distribution function value calculated using a normal
+     *         approximation
+     * @throws MathException if an error occurs computing the normal
+     *             approximation
      */
     public double normalApproximateProbability(int x) throws MathException {
         // calculate the probability using half-correction
@@ -151,9 +157,9 @@
 
     /**
      * Access the domain value lower bound, based on <code>p</code>, used to
-     * bracket a CDF root.  This method is used by
+     * bracket a CDF root. This method is used by
      * {@link #inverseCumulativeProbability(double)} to find critical values.
-     *
+     * 
      * @param p the desired probability for the critical value
      * @return domain lower bound
      */
@@ -164,9 +170,9 @@
 
     /**
      * Access the domain value upper bound, based on <code>p</code>, used to
-     * bracket a CDF root.  This method is used by
+     * bracket a CDF root. This method is used by
      * {@link #inverseCumulativeProbability(double)} to find critical values.
-     *
+     * 
      * @param p the desired probability for the critical value
      * @return domain upper bound
      */
@@ -176,9 +182,10 @@
     }
 
     /**
-     * Modify the normal distribution used to compute normal approximations.
-     * The caller is responsible for insuring the normal distribution has the
-     * proper parameter settings.
+     * Modify the normal distribution used to compute normal approximations. The
+     * caller is responsible for insuring the normal distribution has the proper
+     * parameter settings.
+     * 
      * @param value the new distribution
      * @since 1.2
      */

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/SaddlePointExpansion.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/SaddlePointExpansion.java?rev=830745&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/SaddlePointExpansion.java
(added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/SaddlePointExpansion.java
Wed Oct 28 19:59:21 2009
@@ -0,0 +1,185 @@
+package org.apache.commons.math.distribution;
+
+import org.apache.commons.math.special.Gamma;
+
+/**
+ * <p>
+ * Utility class used by various distributions to accurately compute their
+ * respective probability mass functions. The implementation for this class is
+ * based on the Catherine Loader's <a target="_blank"
+ * href="http://www.herine.net/stat/software/dbinom.html">dbinom</a> routines.
+ * </p>
+ * <p>
+ * This class is not intended to be called directly.
+ * </p>
+ * <p>
+ * References:
+ * <ol>
+ * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
+ * Probabilities.". <a target="_blank"
+ * href="http://www.herine.net/stat/papers/dbinom.pdf">
+ * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
+ * </ol>
+ * </p>
+ * 
+ * @since 1.2
+ * @version $Revision: 1.3 $ $Date: 2007/11/18 23:51:21 $
+ */
+final class SaddlePointExpansion {
+
+    /** 2 &#960;. */
+    private static double PI_2 = 2.0 * Math.PI;
+
+    /** 1/2 * log(2 &#960;). */
+    private static double HALF_LOG_2_PI = 0.5 * Math.log(PI_2);
+
+    /** exact striling expansion error for certain values. */
+    private static final double[] EXACT_STIRLING_ERRORS = { 0.0, /* 0.0 */
+    0.1534264097200273452913848, /* 0.5 */
+    0.0810614667953272582196702, /* 1.0 */
+    0.0548141210519176538961390, /* 1.5 */
+    0.0413406959554092940938221, /* 2.0 */
+    0.03316287351993628748511048, /* 2.5 */
+    0.02767792568499833914878929, /* 3.0 */
+    0.02374616365629749597132920, /* 3.5 */
+    0.02079067210376509311152277, /* 4.0 */
+    0.01848845053267318523077934, /* 4.5 */
+    0.01664469118982119216319487, /* 5.0 */
+    0.01513497322191737887351255, /* 5.5 */
+    0.01387612882307074799874573, /* 6.0 */
+    0.01281046524292022692424986, /* 6.5 */
+    0.01189670994589177009505572, /* 7.0 */
+    0.01110455975820691732662991, /* 7.5 */
+    0.010411265261972096497478567, /* 8.0 */
+    0.009799416126158803298389475, /* 8.5 */
+    0.009255462182712732917728637, /* 9.0 */
+    0.008768700134139385462952823, /* 9.5 */
+    0.008330563433362871256469318, /* 10.0 */
+    0.007934114564314020547248100, /* 10.5 */
+    0.007573675487951840794972024, /* 11.0 */
+    0.007244554301320383179543912, /* 11.5 */
+    0.006942840107209529865664152, /* 12.0 */
+    0.006665247032707682442354394, /* 12.5 */
+    0.006408994188004207068439631, /* 13.0 */
+    0.006171712263039457647532867, /* 13.5 */
+    0.005951370112758847735624416, /* 14.0 */
+    0.005746216513010115682023589, /* 14.5 */
+    0.005554733551962801371038690 /* 15.0 */
+    };
+
+    /**
+     * Default constructor.
+     */
+    private SaddlePointExpansion() {
+        super();
+    }
+
+    /**
+     * Compute the error of Stirling's series at the given value.
+     * <p>
+     * References:
+     * <ol>
+     * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web
+     * Resource. <a target="_blank"
+     * href="http://mathworld.wolfram.com/StirlingsSeries.html">
+     * http://mathworld.wolfram.com/StirlingsSeries.html</a></li>
+     * </ol>
+     * </p>
+     * 
+     * @param z the value.
+     * @return the Striling's series error.
+     */
+    static double getStirlingError(double z) {
+        double ret;
+        if (z < 15.0) {
+            double z2 = 2.0 * z;
+            if (Math.floor(z2) == z2) {
+                ret = EXACT_STIRLING_ERRORS[(int) z2];
+            } else {
+                ret = Gamma.logGamma(z + 1.0) - (z + 0.5) * Math.log(z) + z
+                        - HALF_LOG_2_PI;
+            }
+        } else {
+            double z2 = z * z;
+            ret = (0.083333333333333333333 - (0.00277777777777777777778 - (0.00079365079365079365079365
- (0.000595238095238095238095238 - 0.0008417508417508417508417508 / z2)
+                    / z2)
+                    / z2)
+                    / z2)
+                    / z;
+        }
+        return ret;
+    }
+
+    /**
+     * A part of the deviance portion of the saddle point approximation.
+     * <p>
+     * References:
+     * <ol>
+     * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
+     * Probabilities.". <a target="_blank"
+     * href="http://www.herine.net/stat/papers/dbinom.pdf">
+     * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
+     * </ol>
+     * </p>
+     * 
+     * @param x the x value.
+     * @param mu the average.
+     * @return a part of the deviance.
+     */
+    static double getDeviancePart(double x, double mu) {
+        double ret;
+        if (Math.abs(x - mu) < 0.1 * (x + mu)) {
+            double d = (x - mu);
+            double v = d / (x + mu);
+            double s1 = v * d;
+            double s = Double.NaN;
+            double ej = 2.0 * x * v;
+            v = v * v;
+            int j = 1;
+            while (s1 != s) {
+                s = s1;
+                ej *= v;
+                s1 = s + ej / ((j * 2) + 1);
+                ++j;
+            }
+            ret = s1;
+        } else {
+            ret = x * Math.log(x / mu) + mu - x;
+        }
+        return ret;
+    }
+
+    /**
+     * Compute the PMF for a binomial distribution using the saddle point
+     * expansion.
+     * 
+     * @param x the value at which the probability is evaluated.
+     * @param n the number of trials.
+     * @param p the probability of success.
+     * @param q the probability of failure (1 - p).
+     * @return log(p(x)).
+     */
+    static double logBinomialProbability(int x, int n, double p, double q) {
+        double ret;
+        if (x == 0) {
+            if (p < 0.1) {
+                ret = -getDeviancePart(n, n * q) - n * p;
+            } else {
+                ret = n * Math.log(q);
+            }
+        } else if (x == n) {
+            if (q < 0.1) {
+                ret = -getDeviancePart(n, n * p) - n * q;
+            } else {
+                ret = n * Math.log(p);
+            }
+        } else {
+            ret = getStirlingError(n) - getStirlingError(x)
+                    - getStirlingError(n - x) - getDeviancePart(x, n * p)
+                    - getDeviancePart(n - x, n * q);
+            double f = (PI_2 * x * (n - x)) / n;
+            ret = -0.5 * Math.log(f) + ret;
+        }
+        return ret;
+    }
+}

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=830745&r1=830744&r2=830745&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Wed Oct 28 19:59:21 2009
@@ -39,6 +39,10 @@
   </properties>
   <body>
     <release version="2.1" date="TBD" description="TBD">
+      <action dev="brentworden" type="update" issue="MATH-311" due-to="Nipun Jawalkar">
+        Changed probability calculations for Binomial, Poisson, and Hypergeometric
+        distributions to use Catherine Loader's saddle point approximations.
+      </action>
       <action dev="psteitz" type="fix" issue="MATH-306" due-to="Joerg Huber">
         Removed dead code from Complex#divide.
       </action>

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/distribution/BinomialDistributionTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/distribution/BinomialDistributionTest.java?rev=830745&r1=830744&r2=830745&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/distribution/BinomialDistributionTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/distribution/BinomialDistributionTest.java
Wed Oct 28 19:59:21 2009
@@ -1,57 +1,57 @@
 /*
  * 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.
+ * 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;
 
 /**
- * Test cases for BinomialDistribution.
- * Extends IntegerDistributionAbstractTest.  See class javadoc for
- * IntegerDistributionAbstractTest for details.
- *
- * @version $Revision$ $Date$
+ * Test cases for BinomialDistribution. Extends IntegerDistributionAbstractTest.
+ * See class javadoc for IntegerDistributionAbstractTest for details.
+ * 
+ * @version $Revision$ $Date: 2009-09-05 12:36:48 -0500 (Sat, 05 Sep
+ *          2009) $
  */
 public class BinomialDistributionTest extends IntegerDistributionAbstractTest {
 
     /**
      * Constructor for BinomialDistributionTest.
+     * 
      * @param name
      */
     public BinomialDistributionTest(String name) {
         super(name);
     }
 
-    //-------------- Implementations for abstract methods -----------------------
+    // -------------- Implementations for abstract methods
+    // -----------------------
 
     /** Creates the default discrete distribution instance to use in tests. */
     @Override
     public IntegerDistribution makeDistribution() {
-        return new BinomialDistributionImpl(10,0.70);
+        return new BinomialDistributionImpl(10, 0.70);
     }
 
     /** Creates the default probability density test input values */
     @Override
     public int[] makeDensityTestPoints() {
-        return new int[] {-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
+        return new int[] { -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
     }
 
     /** Creates the default probability density test expected values */
     @Override
     public double[] makeDensityTestValues() {
-        return new double[] {0d, 0.0000d, 0.0001d, 0.0014d, 0.0090d, 0.0368d, 0.1029d,
-                0.2001d, 0.2668d, 0.2335d, 0.1211d, 0.0282d, 0d};
+        return new double[] { 0d, 0.0000059049d, 0.000137781d, 0.0014467d,
+                0.00900169d, 0.0367569d, 0.102919d, 0.200121d, 0.266828d,
+                0.233474d, 0.121061d, 0.0282475d, 0d };
     }
 
     /** Creates the default cumulative probability density test input values */
@@ -63,48 +63,51 @@
     /** Creates the default cumulative probability density test expected values */
     @Override
     public double[] makeCumulativeTestValues() {
-        return new double[] {0d, 0.0000d, 0.0001d, 0.0016d, 0.0106d, 0.0473d,
-                0.1503d, 0.3504d, 0.6172d, 0.8507d, 0.9718d, 1d, 1d};
-        }
+        return new double[] { 0d, 0.0000d, 0.0001d, 0.0016d, 0.0106d, 0.0473d,
+                0.1503d, 0.3504d, 0.6172d, 0.8507d, 0.9718d, 1d, 1d };
+    }
 
     /** Creates the default inverse cumulative probability test input values */
     @Override
     public double[] makeInverseCumulativeTestPoints() {
-        return new double[] {0, 0.001d, 0.010d, 0.025d, 0.050d, 0.100d, 0.999d,
-                0.990d, 0.975d, 0.950d, 0.900d,1};
-        }
+        return new double[] { 0, 0.001d, 0.010d, 0.025d, 0.050d, 0.100d,
+                0.999d, 0.990d, 0.975d, 0.950d, 0.900d, 1 };
+    }
 
-    /** Creates the default inverse cumulative probability density test expected values */
+    /**
+     * Creates the default inverse cumulative probability density test expected
+     * values
+     */
     @Override
     public int[] makeInverseCumulativeTestValues() {
-        return new int[] {-1, 1, 2, 3, 4, 4, 9, 9, 9, 8, 8, Integer.MAX_VALUE};
+        return new int[] { -1, 1, 2, 3, 4, 4, 9, 9, 9, 8, 8, Integer.MAX_VALUE };
     }
 
-    //----------------- Additional test cases ---------------------------------
+    // ----------------- Additional test cases ---------------------------------
 
-    /** Test degenerate case p = 0   */
+    /** Test degenerate case p = 0 */
     public void testDegenerate0() throws Exception {
-        setDistribution(new BinomialDistributionImpl(5,0.0d));
-        setCumulativeTestPoints(new int[] {-1, 0, 1, 5, 10 });
-        setCumulativeTestValues(new double[] {0d, 1d, 1d, 1d, 1d});
-        setDensityTestPoints(new int[] {-1, 0, 1, 10, 11});
-        setDensityTestValues(new double[] {0d, 1d, 0d, 0d, 0d});
-        setInverseCumulativeTestPoints(new double[] {0.1d, 0.5d});
-        setInverseCumulativeTestValues(new int[] {-1, -1});
+        setDistribution(new BinomialDistributionImpl(5, 0.0d));
+        setCumulativeTestPoints(new int[] { -1, 0, 1, 5, 10 });
+        setCumulativeTestValues(new double[] { 0d, 1d, 1d, 1d, 1d });
+        setDensityTestPoints(new int[] { -1, 0, 1, 10, 11 });
+        setDensityTestValues(new double[] { 0d, 1d, 0d, 0d, 0d });
+        setInverseCumulativeTestPoints(new double[] { 0.1d, 0.5d });
+        setInverseCumulativeTestValues(new int[] { -1, -1 });
         verifyDensities();
         verifyCumulativeProbabilities();
         verifyInverseCumulativeProbabilities();
     }
 
-    /** Test degenerate case p = 1   */
+    /** Test degenerate case p = 1 */
     public void testDegenerate1() throws Exception {
-        setDistribution(new BinomialDistributionImpl(5,1.0d));
-        setCumulativeTestPoints(new int[] {-1, 0, 1, 2, 5, 10 });
-        setCumulativeTestValues(new double[] {0d, 0d, 0d, 0d, 1d, 1d});
-        setDensityTestPoints(new int[] {-1, 0, 1, 2, 5, 10});
-        setDensityTestValues(new double[] {0d, 0d, 0d, 0d, 1d, 0d});
-        setInverseCumulativeTestPoints(new double[] {0.1d, 0.5d});
-        setInverseCumulativeTestValues(new int[] {4, 4});
+        setDistribution(new BinomialDistributionImpl(5, 1.0d));
+        setCumulativeTestPoints(new int[] { -1, 0, 1, 2, 5, 10 });
+        setCumulativeTestValues(new double[] { 0d, 0d, 0d, 0d, 1d, 1d });
+        setDensityTestPoints(new int[] { -1, 0, 1, 2, 5, 10 });
+        setDensityTestValues(new double[] { 0d, 0d, 0d, 0d, 1d, 0d });
+        setInverseCumulativeTestPoints(new double[] { 0.1d, 0.5d });
+        setInverseCumulativeTestValues(new int[] { 4, 4 });
         verifyDensities();
         verifyCumulativeProbabilities();
         verifyInverseCumulativeProbabilities();



Mime
View raw message