View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.distribution.discrete;
23  
24  import org.hipparchus.special.Gamma;
25  import org.hipparchus.util.FastMath;
26  import org.hipparchus.util.MathUtils;
27  import org.hipparchus.util.Precision;
28  
29  /**
30   * Utility class used by various distributions to accurately compute their
31   * respective probability mass functions. The implementation for this class is
32   * based on the Catherine Loader's <a target="_blank"
33   * href="http://www.herine.net/stat/software/dbinom.html">dbinom</a> routines.
34   * <p>
35   * This class is not intended to be called directly.
36   * <p>
37   * References:
38   * <ol>
39   * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
40   * Probabilities.". <a target="_blank"
41   * href="http://www.herine.net/stat/papers/dbinom.pdf">
42   * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
43   * </ol>
44   */
45  final class SaddlePointExpansion {
46  
47      /** 1/2 * log(2 &#960;). */
48      private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(MathUtils.TWO_PI);
49  
50      /** exact Stirling expansion error for certain values. */
51      private static final double[] EXACT_STIRLING_ERRORS = {
52          0.0,                           /* 0.0 */
53          0.1534264097200273452913848,   /* 0.5 */
54          0.0810614667953272582196702,   /* 1.0 */
55          0.0548141210519176538961390,   /* 1.5 */
56          0.0413406959554092940938221,   /* 2.0 */
57          0.03316287351993628748511048,  /* 2.5 */
58          0.02767792568499833914878929,  /* 3.0 */
59          0.02374616365629749597132920,  /* 3.5 */
60          0.02079067210376509311152277,  /* 4.0 */
61          0.01848845053267318523077934,  /* 4.5 */
62          0.01664469118982119216319487,  /* 5.0 */
63          0.01513497322191737887351255,  /* 5.5 */
64          0.01387612882307074799874573,  /* 6.0 */
65          0.01281046524292022692424986,  /* 6.5 */
66          0.01189670994589177009505572,  /* 7.0 */
67          0.01110455975820691732662991,  /* 7.5 */
68          0.010411265261972096497478567, /* 8.0 */
69          0.009799416126158803298389475, /* 8.5 */
70          0.009255462182712732917728637, /* 9.0 */
71          0.008768700134139385462952823, /* 9.5 */
72          0.008330563433362871256469318, /* 10.0 */
73          0.007934114564314020547248100, /* 10.5 */
74          0.007573675487951840794972024, /* 11.0 */
75          0.007244554301320383179543912, /* 11.5 */
76          0.006942840107209529865664152, /* 12.0 */
77          0.006665247032707682442354394, /* 12.5 */
78          0.006408994188004207068439631, /* 13.0 */
79          0.006171712263039457647532867, /* 13.5 */
80          0.005951370112758847735624416, /* 14.0 */
81          0.005746216513010115682023589, /* 14.5 */
82          0.005554733551962801371038690  /* 15.0 */
83      };
84  
85      /**
86       * Default constructor.
87       */
88      private SaddlePointExpansion() {}
89  
90      /**
91       * Compute the error of Stirling's series at the given value.
92       * <p>
93       * References:
94       * <ol>
95       * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web
96       * Resource. <a target="_blank"
97       * href="http://mathworld.wolfram.com/StirlingsSeries.html">
98       * http://mathworld.wolfram.com/StirlingsSeries.html</a></li>
99       * </ol>
100      *
101      * @param z the value.
102      * @return the Striling's series error.
103      */
104     static double getStirlingError(double z) {
105         if (z < 15.0) {
106             double z2 = 2.0 * z;
107             if (Precision.isMathematicalInteger(z2)) {
108                 return EXACT_STIRLING_ERRORS[(int) z2];
109             } else {
110                 return Gamma.logGamma(z + 1.0) - (z + 0.5) * FastMath.log(z) +
111                        z - HALF_LOG_2_PI;
112             }
113         } else {
114             double z2 = z * z;
115             return (0.083333333333333333333 -
116                            (0.00277777777777777777778 -
117                              (0.00079365079365079365079365 -
118                                (0.000595238095238095238095238 -
119                                  0.0008417508417508417508417508 /
120                                    z2) / z2) / z2) / z2) / z;
121         }
122     }
123 
124     /**
125      * A part of the deviance portion of the saddle point approximation.
126      * <p>
127      * References:
128      * <ol>
129      * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
130      * Probabilities.". <a target="_blank"
131      * href="http://www.herine.net/stat/papers/dbinom.pdf">
132      * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
133      * </ol>
134      *
135      * @param x the x value.
136      * @param mu the average.
137      * @return a part of the deviance.
138      */
139     static double getDeviancePart(double x, double mu) {
140         if (FastMath.abs(x - mu) < 0.1 * (x + mu)) {
141             double d = x - mu;
142             double v = d / (x + mu);
143             double s1 = v * d;
144             double s = Double.NaN;
145             double ej = 2.0 * x * v;
146             v *= v;
147             int j = 1;
148             while (s1 != s) {
149                 s = s1;
150                 ej *= v;
151                 s1 = s + ej / ((j * 2) + 1);
152                 ++j;
153             }
154             return s1;
155         } else {
156             return x * FastMath.log(x / mu) + mu - x;
157         }
158     }
159 
160     /**
161      * Compute the logarithm of the PMF for a binomial distribution
162      * using the saddle point expansion.
163      *
164      * @param x the value at which the probability is evaluated.
165      * @param n the number of trials.
166      * @param p the probability of success.
167      * @param q the probability of failure (1 - p).
168      * @return log(p(x)).
169      */
170     static double logBinomialProbability(int x, int n, double p, double q) {
171         if (n == 0) {
172             return x == 0 ? 0d : Double.NEGATIVE_INFINITY;
173         }
174         if (x == 0) {
175             if (p < 0.1) {
176                 return -getDeviancePart(n, n * q) - n * p;
177             } else {
178                 return n * FastMath.log(q);
179             }
180         } else if (x == n) {
181             if (q < 0.1) {
182                 return -getDeviancePart(n, n * p) - n * q;
183             } else {
184                 return n * FastMath.log(p);
185             }
186         } else {
187             double ret = getStirlingError(n) - getStirlingError(x) -
188                          getStirlingError(n - x) - getDeviancePart(x, n * p) -
189                          getDeviancePart(n - x, n * q);
190             double f = (MathUtils.TWO_PI * x * (n - x)) / n;
191             ret = -0.5 * FastMath.log(f) + ret;
192             return ret;
193         }
194     }
195 }