source: src/main/java/agents/anac/y2019/harddealer/math3/distribution/SaddlePointExpansion.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents

  • Property svn:executable set to *
File size: 6.9 KB
Line 
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 */
17package agents.anac.y2019.harddealer.math3.distribution;
18
19import agents.anac.y2019.harddealer.math3.special.Gamma;
20import agents.anac.y2019.harddealer.math3.util.FastMath;
21import 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 */
45final 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 = { 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}
Note: See TracBrowser for help on using the repository browser.