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 agents.anac.y2019.harddealer.math3.distribution;
|
---|
18 |
|
---|
19 | import agents.anac.y2019.harddealer.math3.special.Gamma;
|
---|
20 | import agents.anac.y2019.harddealer.math3.util.FastMath;
|
---|
21 | import agents.anac.y2019.harddealer.math3.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 | */
|
---|
45 | final class SaddlePointExpansion {
|
---|
46 |
|
---|
47 | /** 1/2 * log(2 π). */
|
---|
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 = { 0.0, /* 0.0 */
|
---|
52 | 0.1534264097200273452913848, /* 0.5 */
|
---|
53 | 0.0810614667953272582196702, /* 1.0 */
|
---|
54 | 0.0548141210519176538961390, /* 1.5 */
|
---|
55 | 0.0413406959554092940938221, /* 2.0 */
|
---|
56 | 0.03316287351993628748511048, /* 2.5 */
|
---|
57 | 0.02767792568499833914878929, /* 3.0 */
|
---|
58 | 0.02374616365629749597132920, /* 3.5 */
|
---|
59 | 0.02079067210376509311152277, /* 4.0 */
|
---|
60 | 0.01848845053267318523077934, /* 4.5 */
|
---|
61 | 0.01664469118982119216319487, /* 5.0 */
|
---|
62 | 0.01513497322191737887351255, /* 5.5 */
|
---|
63 | 0.01387612882307074799874573, /* 6.0 */
|
---|
64 | 0.01281046524292022692424986, /* 6.5 */
|
---|
65 | 0.01189670994589177009505572, /* 7.0 */
|
---|
66 | 0.01110455975820691732662991, /* 7.5 */
|
---|
67 | 0.010411265261972096497478567, /* 8.0 */
|
---|
68 | 0.009799416126158803298389475, /* 8.5 */
|
---|
69 | 0.009255462182712732917728637, /* 9.0 */
|
---|
70 | 0.008768700134139385462952823, /* 9.5 */
|
---|
71 | 0.008330563433362871256469318, /* 10.0 */
|
---|
72 | 0.007934114564314020547248100, /* 10.5 */
|
---|
73 | 0.007573675487951840794972024, /* 11.0 */
|
---|
74 | 0.007244554301320383179543912, /* 11.5 */
|
---|
75 | 0.006942840107209529865664152, /* 12.0 */
|
---|
76 | 0.006665247032707682442354394, /* 12.5 */
|
---|
77 | 0.006408994188004207068439631, /* 13.0 */
|
---|
78 | 0.006171712263039457647532867, /* 13.5 */
|
---|
79 | 0.005951370112758847735624416, /* 14.0 */
|
---|
80 | 0.005746216513010115682023589, /* 14.5 */
|
---|
81 | 0.005554733551962801371038690 /* 15.0 */
|
---|
82 | };
|
---|
83 |
|
---|
84 | /**
|
---|
85 | * Default constructor.
|
---|
86 | */
|
---|
87 | private SaddlePointExpansion() {
|
---|
88 | super();
|
---|
89 | }
|
---|
90 |
|
---|
91 | /**
|
---|
92 | * Compute the error of Stirling's series at the given value.
|
---|
93 | * <p>
|
---|
94 | * References:
|
---|
95 | * <ol>
|
---|
96 | * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web
|
---|
97 | * Resource. <a target="_blank"
|
---|
98 | * href="http://mathworld.wolfram.com/StirlingsSeries.html">
|
---|
99 | * http://mathworld.wolfram.com/StirlingsSeries.html</a></li>
|
---|
100 | * </ol>
|
---|
101 | * </p>
|
---|
102 | *
|
---|
103 | * @param z the value.
|
---|
104 | * @return the Striling's series error.
|
---|
105 | */
|
---|
106 | static double getStirlingError(double z) {
|
---|
107 | double ret;
|
---|
108 | if (z < 15.0) {
|
---|
109 | double z2 = 2.0 * z;
|
---|
110 | if (FastMath.floor(z2) == z2) {
|
---|
111 | ret = EXACT_STIRLING_ERRORS[(int) z2];
|
---|
112 | } else {
|
---|
113 | ret = Gamma.logGamma(z + 1.0) - (z + 0.5) * FastMath.log(z) +
|
---|
114 | z - HALF_LOG_2_PI;
|
---|
115 | }
|
---|
116 | } else {
|
---|
117 | double z2 = z * z;
|
---|
118 | ret = (0.083333333333333333333 -
|
---|
119 | (0.00277777777777777777778 -
|
---|
120 | (0.00079365079365079365079365 -
|
---|
121 | (0.000595238095238095238095238 -
|
---|
122 | 0.0008417508417508417508417508 /
|
---|
123 | z2) / z2) / z2) / z2) / z;
|
---|
124 | }
|
---|
125 | return ret;
|
---|
126 | }
|
---|
127 |
|
---|
128 | /**
|
---|
129 | * A part of the deviance portion of the saddle point approximation.
|
---|
130 | * <p>
|
---|
131 | * References:
|
---|
132 | * <ol>
|
---|
133 | * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
|
---|
134 | * Probabilities.". <a target="_blank"
|
---|
135 | * href="http://www.herine.net/stat/papers/dbinom.pdf">
|
---|
136 | * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
|
---|
137 | * </ol>
|
---|
138 | * </p>
|
---|
139 | *
|
---|
140 | * @param x the x value.
|
---|
141 | * @param mu the average.
|
---|
142 | * @return a part of the deviance.
|
---|
143 | */
|
---|
144 | static double getDeviancePart(double x, double mu) {
|
---|
145 | double ret;
|
---|
146 | if (FastMath.abs(x - mu) < 0.1 * (x + mu)) {
|
---|
147 | double d = x - mu;
|
---|
148 | double v = d / (x + mu);
|
---|
149 | double s1 = v * d;
|
---|
150 | double s = Double.NaN;
|
---|
151 | double ej = 2.0 * x * v;
|
---|
152 | v *= v;
|
---|
153 | int j = 1;
|
---|
154 | while (s1 != s) {
|
---|
155 | s = s1;
|
---|
156 | ej *= v;
|
---|
157 | s1 = s + ej / ((j * 2) + 1);
|
---|
158 | ++j;
|
---|
159 | }
|
---|
160 | ret = s1;
|
---|
161 | } else {
|
---|
162 | ret = x * FastMath.log(x / mu) + mu - x;
|
---|
163 | }
|
---|
164 | return ret;
|
---|
165 | }
|
---|
166 |
|
---|
167 | /**
|
---|
168 | * Compute the logarithm of the PMF for a binomial distribution
|
---|
169 | * using the saddle point expansion.
|
---|
170 | *
|
---|
171 | * @param x the value at which the probability is evaluated.
|
---|
172 | * @param n the number of trials.
|
---|
173 | * @param p the probability of success.
|
---|
174 | * @param q the probability of failure (1 - p).
|
---|
175 | * @return log(p(x)).
|
---|
176 | */
|
---|
177 | static double logBinomialProbability(int x, int n, double p, double q) {
|
---|
178 | double ret;
|
---|
179 | if (x == 0) {
|
---|
180 | if (p < 0.1) {
|
---|
181 | ret = -getDeviancePart(n, n * q) - n * p;
|
---|
182 | } else {
|
---|
183 | ret = n * FastMath.log(q);
|
---|
184 | }
|
---|
185 | } else if (x == n) {
|
---|
186 | if (q < 0.1) {
|
---|
187 | ret = -getDeviancePart(n, n * p) - n * q;
|
---|
188 | } else {
|
---|
189 | ret = n * FastMath.log(p);
|
---|
190 | }
|
---|
191 | } else {
|
---|
192 | ret = getStirlingError(n) - getStirlingError(x) -
|
---|
193 | getStirlingError(n - x) - getDeviancePart(x, n * p) -
|
---|
194 | getDeviancePart(n - x, n * q);
|
---|
195 | double f = (MathUtils.TWO_PI * x * (n - x)) / n;
|
---|
196 | ret = -0.5 * FastMath.log(f) + ret;
|
---|
197 | }
|
---|
198 | return ret;
|
---|
199 | }
|
---|
200 | }
|
---|