wfmath 1.0.3
A math library for the Worldforge system.
probability.cpp
1// probability.cpp (probability and statistics implementation)
2//
3// The WorldForge Project
4// Copyright (C) 2002 The WorldForge Project
5//
6// This program is free software; you can redistribute it and/or modify
7// it under the terms of the GNU General Public License as published by
8// the Free Software Foundation; either version 2 of the License, or
9// (at your option) any later version.
10//
11// This program is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14// GNU General Public License for more details.
15//
16// You should have received a copy of the GNU General Public License
17// along with this program; if not, write to the Free Software
18// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
19//
20// For information about WorldForge and its authors, please contact
21// the Worldforge Web Site at http://www.worldforge.org.
22
23// Author: Ron Steinke
24// Created: 2002-1-23
25
26#include <wfmath/probability.h>
27
28#include <wfmath/const.h>
29
30#include <cmath>
31
32#include <cassert>
33
34namespace WFMath {
35
36template<typename FloatT>
37static FloatT LogPoisson(FloatT mean, unsigned int step);
38template<typename FloatT>
39static FloatT IncompleteGamma(FloatT a, FloatT z);
40template<typename FloatT>
41static FloatT IncompleteGammaNoPrefactor(FloatT a, FloatT z);
42template<typename FloatT>
43static FloatT IncompleteGammaComplement(FloatT a, FloatT z);
44template<typename FloatT>
45static FloatT IncompleteGammaComplementNoPrefactor(FloatT a, FloatT z);
46
47// Making this an int makes LogFactorial faster
48static const unsigned int GammaCutoff = 10;
49
50template<typename FloatT>
51FloatT GaussianConditional(FloatT mean, FloatT stddev, FloatT val)
52{
53 assert(stddev != 0);
54
55 FloatT diff = val - mean;
56 FloatT diffnorm = diff / stddev;
57 FloatT diffsqr_over_two = diffnorm * diffnorm / 2;
58
59 /* Make sure round off error in Sqrt3 doesn't hit
60 * assert() in IncompleteGammaComplementNoPrefactor()
61 */
62 static const FloatT half = 0.5;
63 if(diffnorm < numeric_constants<FloatT>::sqrt3() +
64 10 * std::numeric_limits<FloatT>::epsilon()) {
65 FloatT erfc_norm = IncompleteGammaComplement(half, diffsqr_over_two);
66
67 FloatT normalization = (diffnorm > 0) ? (erfc_norm / 2) : (1 - erfc_norm / 2);
68
69 return Gaussian(mean, stddev, val) / normalization;
70 }
71 static const FloatT two = 2.0;
72
73 return two / (std::fabs(diff)
74 * IncompleteGammaComplementNoPrefactor<FloatT>(half, diffsqr_over_two));
75}
76
77template<typename FloatT>
78FloatT Gaussian(FloatT mean, FloatT stddev, FloatT val)
79{
80 assert(stddev != 0);
81
82 FloatT diff = (mean - val) / stddev;
83
84 return std::exp(-(diff * diff) / 2) / (std::fabs(stddev) *
87}
88
89template<typename FloatT>
90FloatT PoissonConditional(FloatT mean, unsigned int step)
91{
92 assert(mean >= 0);
93
94 if(mean == 0) // Funky limit, but allow it
95 return (step == 0) ? 1 : 0;
96
97 if(step == 0)
98 return std::exp(-mean);
99
100 if(mean > step + 1)
101 return Poisson(mean, step) /
102 IncompleteGamma(static_cast<FloatT>(step), mean);
103
104 static const FloatT one = 1.0;
105 return one / IncompleteGammaNoPrefactor(static_cast<FloatT>(step), mean);
106}
107
108template<typename FloatT>
109FloatT Poisson(FloatT mean, unsigned int step)
110{
111 assert(mean >= 0);
112
113 if(mean == 0) // Funky limit, but allow it
114 return (step == 0) ? 1 : 0;
115
116 return std::exp(LogPoisson(mean, step));
117}
118
119template<typename FloatT>
120static FloatT LogPoisson(FloatT mean, unsigned int step)
121{
122 assert(mean > 0);
123
124 if(step == 0)
125 return -mean;
126
127 FloatT first = static_cast<FloatT>(step) * std::log(mean);
128 FloatT second = mean + LogFactorial<FloatT>(step);
129
130 assert("LogFactorial() always returns positive" && second > 0);
131
132 FloatT ans = first - second;
133
134 // can only get probability == 1 for step == mean == 0
135 assert("must get probability < 1" && ans < 0);
136
137 return ans;
138}
139
140template<typename FloatT>
141FloatT Factorial(unsigned int n)
142{
143 if(n == 0 || n == 1)
144 return 1;
145
146 if(n < GammaCutoff) {
147 FloatT ans = static_cast<FloatT>(n);
148 while(--n > 1) // Don't need to multiply by 1
149 ans *= static_cast<FloatT>(n);
150 return ans;
151 }
152 else
153 return std::exp(LogGamma(static_cast<FloatT>(n + 1)));
154}
155
156template<typename FloatT>
157FloatT LogFactorial(unsigned int n)
158{
159 if(n == 0 || n == 1)
160 return 0; // ln(0!) = ln(1!) = ln(1) = 0
161
162 if(n < GammaCutoff) {
163 FloatT ans = static_cast<FloatT>(n);
164 while(--n > 1) // Don't need to multiply by 1
165 ans *= static_cast<FloatT>(n);
166 return std::log(ans);
167 }
168 else
169 return LogGamma(static_cast<FloatT>(n + 1));
170}
171
172template<typename FloatT>
173FloatT Gamma(FloatT z)
174{
175 if(z >= 0.5)
176 return std::exp(LogGamma(z));
177 else
179 std::exp(-LogGamma(1 - z)) /
180 std::sin(numeric_constants<FloatT>::pi() * z);
181}
182
183template<typename FloatT>
184FloatT LogGamma(FloatT z)
185{
186 /* This will return nan or something when z is a non-positive
187 * integer (from trying to take log(0)), but that's what
188 * should happen anyway.
189 */
190 static const FloatT one = 1.0;
191
192 if(z < 0.5)
194 LogGamma<FloatT>(one - z) -
195 std::log(std::fabs(std::sin(numeric_constants<FloatT>::pi() * z)));
196
197 if(z == 0.5) // special case for Gaussian
199
200 if(z == 1 || z == 2) // 0! and 1!
201 return 0;
202
203 FloatT log_shift;
204
205 if(z < GammaCutoff) {
206 FloatT shift = 1;
207 while(z < GammaCutoff)
208 shift *= z++;
209 log_shift = std::log(std::fabs(shift));
210 }
211 else
212 log_shift = 0;
213
214 // Stirling approximation (see Gradshteyn + Ryzhik, Table of Integrals,
215 // Series, and Products, fifth edition, formula 8.344 for a specific formula)
216
217 static const FloatT half = 0.5, two = 2.0;
218
219 FloatT ans = (z - half) * std::log(z) - log_shift - z +
222
223 // coeffs[i] is the 2*(i+1)th Bernoulli number, divided by (2i + 1)*(2i + 2)
224 static const FloatT coeffs[] = {1.0/12.0, -1.0/360.0,
225 1.0/1260.0, -1.0/1620.0, 5.0/5940.0,
226 -691.0/360360.0, 7.0/1092.0, -3617.0/122400.0,
227 43867.0/244188.0, -174661.0/125400.0, 854513.0/63756.0,};
228 static const int num_coeffs = sizeof(coeffs) / sizeof(FloatT);
229
230 FloatT z_power = 1/z;
231 FloatT z_to_minus_two = z_power * z_power;
232 FloatT small_enough = std::fabs(ans) * std::numeric_limits<FloatT>::epsilon();
233 int i;
234
235 for(i = 0; i < num_coeffs; ++i) {
236 FloatT next_term = coeffs[i] * z_power;
237 ans += next_term;
238 if(std::fabs(next_term) < small_enough)
239 break;
240 z_power *= z_to_minus_two;
241 }
242
243 assert(i < num_coeffs); // If someone hits this, tell me and I'll add more terms,
244 // worst case is for n = cutoff = 10, which should work
245 // for DBL_EPSILON > 1.048e-21
246
247 return ans;
248}
249
250template<typename FloatT>
251static FloatT IncompleteGamma(FloatT a, FloatT z)
252{
253 assert(z >= 0 && a >= 0); // We only use this part of the parameter space
254
255 if(a == 0)
256 return 1;
257 else if(z == 0)
258 return 0;
259
260 if(z > a + 1)
261 return 1 - IncompleteGammaComplement(a, z);
262
263 FloatT prefactor = std::exp(a * (std::log(z) + 1) - z - LogGamma(a));
264
265 return IncompleteGammaNoPrefactor(a, z) * prefactor;
266}
267
268template<typename FloatT>
269static FloatT IncompleteGammaNoPrefactor(FloatT a, FloatT z)
270{
271 assert(z <= a + 1 && z > 0 && a > 0);
272
273 // Power series
274
275 FloatT term = 1;
276 FloatT dividend = a;
277
278 FloatT ans = term;
279 while(std::fabs(term / ans) > std::numeric_limits<FloatT>::epsilon()) {
280 term *= z / ++dividend;
281 ans += term;
282 }
283
284 return ans;
285}
286
287template<typename FloatT>
288static FloatT IncompleteGammaComplement(FloatT a, FloatT z)
289{
290 assert(z >= 0 && a >= 0); // We only use this part of the parameter space
291
292 if(a == 0)
293 return 0;
294 else if(z == 0)
295 return 1;
296
297 if(z < a + 1)
298 return 1 - IncompleteGamma(a, z);
299
300 FloatT prefactor = std::exp(a * std::log(z) - z - LogGamma(a));
301
302 return IncompleteGammaComplementNoPrefactor(a, z) * prefactor;
303}
304
305template<typename FloatT>
306static FloatT IncompleteGammaComplementNoPrefactor(FloatT a, FloatT z)
307{
308 assert(z >= a + 1 && a > 0);
309
310 // Continued fraction
311
312 static const FloatT fudge = 1000;
313
314 FloatT b_contrib = z + 1 - a;
315 FloatT a_last, b_last, a_next, b_next;
316 FloatT term = 1;
317 bool last_zero,
318 next_zero = (std::fabs(b_contrib) <= std::numeric_limits<FloatT>::min()
319 * fudge);
320
321 if(next_zero) {
322 a_last = 0;
323 b_last = 1;
324 a_next = 1;
325 b_next = b_contrib;
326 }
327 else {
328 a_last = 0;
329 b_last = 1 / b_contrib;
330 a_next = b_last;
331 b_next = 1;
332 }
333
334 while(true) {
335 FloatT a_tmp = a_next, b_tmp = b_next;
336
337 FloatT a_contrib = term * (a - term);
338 ++term;
339 b_contrib += 2;
340
341 a_next = b_contrib * a_tmp + a_contrib * a_last;
342 b_next = b_contrib * b_tmp + a_contrib * b_last;
343
344 a_last = a_tmp;
345 b_last = b_tmp;
346
347 last_zero = next_zero;
348 next_zero = (std::fabs(b_next) <=
349 std::fabs(a_next) * (std::numeric_limits<FloatT>::min() *
350 fudge));
351
352 if(next_zero)
353 continue; // b_next is about zero
354
355 a_next /= b_next;
356
357 if(!last_zero &&
358 std::fabs(a_next - a_last) < std::fabs(a_last) *
359 std::numeric_limits<FloatT>::epsilon())
360 return a_next; // Can't compare if b_last was zero
361
362 a_last /= b_next;
363 b_last /= b_next;
364 b_next = 1;
365 }
366}
367
368template
369float GaussianConditional<float>(float mean, float stddev, float val);
370template
371float Gaussian<float>(float mean, float stddev, float val);
372template
373float PoissonConditional<float>(float mean, unsigned int step);
374template
375float Poisson<float>(float mean, unsigned int step);
376template
377float LogFactorial<float>(unsigned int n);
378template
379float Factorial<float>(unsigned int n);
380template
381float LogGamma<float>(float z);
382template
383float Gamma<float>(float z);
384
385template
386double GaussianConditional<double>(double mean, double stddev, double val);
387template
388double Gaussian<double>(double mean, double stddev, double val);
389template
390double PoissonConditional<double>(double mean, unsigned int step);
391template
392double Poisson<double>(double mean, unsigned int step);
393template
394double LogFactorial<double>(unsigned int n);
395template
396double Factorial<double>(unsigned int n);
397template
398double LogGamma<double>(double z);
399template
400double Gamma<double>(double z);
401
402} // namespace WFMath
Generic library namespace.
Definition: shape.h:41
FloatT GaussianConditional(FloatT mean, FloatT stddev, FloatT val)
Gives the conditional probability of the Gaussian distribution at position val.
Definition: probability.cpp:51
FloatT Poisson(FloatT mean, unsigned int step)
Gives the value of the Poisson distribution at position step.
FloatT LogGamma(FloatT z)
The natural log of Euler's Gamma function.
FloatT Factorial(unsigned int n)
Gives n!
FloatT LogFactorial(unsigned int n)
Gives the natural log of n!
FloatT Gaussian(FloatT mean, FloatT stddev, FloatT val)
Gives the value of the Gaussian distribution at position val.
Definition: probability.cpp:78
FloatT Gamma(FloatT z)
Euler's Gamma function.
FloatT PoissonConditional(FloatT mean, unsigned int step)
Gives the conditional probability of the Poisson distribution at position step.
Definition: probability.cpp:90
static FloatType pi()
The constant pi.
static FloatType log2()
The natural logarithm of 2.
static FloatType log_pi()
The natural logarithm of pi.
static FloatType sqrt2()
The square root of 2.
static FloatType sqrt_pi()
The square root of pi.