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 * http://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 package org.apache.commons.math.distribution; 18 19 import org.apache.commons.math.special.Gamma; 20 import org.apache.commons.math.util.FastMath; 21 import org.apache.commons.math.util.MathUtils; 22 23 /** 24 * <p> 25 * Utility class used by various distributions to accurately compute their 26 * respective probability mass functions. The implementation for this class is 27 * based on the Catherine Loader's <a target="_blank" 28 * href="http://www.herine.net/stat/software/dbinom.html">dbinom</a> routines. 29 * </p> 30 * <p> 31 * This class is not intended to be called directly. 32 * </p> 33 * <p> 34 * References: 35 * <ol> 36 * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial 37 * Probabilities.". <a target="_blank" 38 * href="http://www.herine.net/stat/papers/dbinom.pdf"> 39 * http://www.herine.net/stat/papers/dbinom.pdf</a></li> 40 * </ol> 41 * </p> 42 * 43 * @since 2.1 44 * @version $Revision$ $Date$ 45 */ 46 final class SaddlePointExpansion { 47 48 /** 1/2 * log(2 π). */ 49 private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(MathUtils.TWO_PI); 50 51 /** exact Stirling expansion error for certain values. */ 52 private static final double[] EXACT_STIRLING_ERRORS = { 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 super(); 90 } 91 92 /** 93 * Compute the error of Stirling's series at the given value. 94 * <p> 95 * References: 96 * <ol> 97 * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web 98 * Resource. <a target="_blank" 99 * href="http://mathworld.wolfram.com/StirlingsSeries.html"> 100 * http://mathworld.wolfram.com/StirlingsSeries.html</a></li> 101 * </ol> 102 * </p> 103 * 104 * @param z the value. 105 * @return the Striling's series error. 106 */ 107 static double getStirlingError(double z) { 108 double ret; 109 if (z < 15.0) { 110 double z2 = 2.0 * z; 111 if (FastMath.floor(z2) == z2) { 112 ret = EXACT_STIRLING_ERRORS[(int) z2]; 113 } else { 114 ret = Gamma.logGamma(z + 1.0) - (z + 0.5) * FastMath.log(z) + 115 z - HALF_LOG_2_PI; 116 } 117 } else { 118 double z2 = z * z; 119 ret = (0.083333333333333333333 - 120 (0.00277777777777777777778 - 121 (0.00079365079365079365079365 - 122 (0.000595238095238095238095238 - 123 0.0008417508417508417508417508 / 124 z2) / z2) / z2) / z2) / z; 125 } 126 return ret; 127 } 128 129 /** 130 * A part of the deviance portion of the saddle point approximation. 131 * <p> 132 * References: 133 * <ol> 134 * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial 135 * Probabilities.". <a target="_blank" 136 * href="http://www.herine.net/stat/papers/dbinom.pdf"> 137 * http://www.herine.net/stat/papers/dbinom.pdf</a></li> 138 * </ol> 139 * </p> 140 * 141 * @param x the x value. 142 * @param mu the average. 143 * @return a part of the deviance. 144 */ 145 static double getDeviancePart(double x, double mu) { 146 double ret; 147 if (FastMath.abs(x - mu) < 0.1 * (x + mu)) { 148 double d = x - mu; 149 double v = d / (x + mu); 150 double s1 = v * d; 151 double s = Double.NaN; 152 double ej = 2.0 * x * v; 153 v = v * v; 154 int j = 1; 155 while (s1 != s) { 156 s = s1; 157 ej *= v; 158 s1 = s + ej / ((j * 2) + 1); 159 ++j; 160 } 161 ret = s1; 162 } else { 163 ret = x * FastMath.log(x / mu) + mu - x; 164 } 165 return ret; 166 } 167 168 /** 169 * Compute the PMF for a binomial distribution using the saddle point 170 * expansion. 171 * 172 * @param x the value at which the probability is evaluated. 173 * @param n the number of trials. 174 * @param p the probability of success. 175 * @param q the probability of failure (1 - p). 176 * @return log(p(x)). 177 */ 178 static double logBinomialProbability(int x, int n, double p, double q) { 179 double ret; 180 if (x == 0) { 181 if (p < 0.1) { 182 ret = -getDeviancePart(n, n * q) - n * p; 183 } else { 184 ret = n * FastMath.log(q); 185 } 186 } else if (x == n) { 187 if (q < 0.1) { 188 ret = -getDeviancePart(n, n * p) - n * q; 189 } else { 190 ret = n * FastMath.log(p); 191 } 192 } else { 193 ret = getStirlingError(n) - getStirlingError(x) - 194 getStirlingError(n - x) - getDeviancePart(x, n * p) - 195 getDeviancePart(n - x, n * q); 196 double f = (MathUtils.TWO_PI * x * (n - x)) / n; 197 ret = -0.5 * FastMath.log(f) + ret; 198 } 199 return ret; 200 } 201 } 202
