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 
34 namespace WFMath {
35 
36 template<typename FloatT>
37 static FloatT LogPoisson(FloatT mean, unsigned int step);
38 template<typename FloatT>
39 static FloatT IncompleteGamma(FloatT a, FloatT z);
40 template<typename FloatT>
41 static FloatT IncompleteGammaNoPrefactor(FloatT a, FloatT z);
42 template<typename FloatT>
43 static FloatT IncompleteGammaComplement(FloatT a, FloatT z);
44 template<typename FloatT>
45 static FloatT IncompleteGammaComplementNoPrefactor(FloatT a, FloatT z);
46 
47 // Making this an int makes LogFactorial faster
48 static const unsigned int GammaCutoff = 10;
49 
50 template<typename FloatT>
51 FloatT 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 
77 template<typename FloatT>
78 FloatT 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 
89 template<typename FloatT>
90 FloatT 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 
108 template<typename FloatT>
109 FloatT 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 
119 template<typename FloatT>
120 static 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 
140 template<typename FloatT>
141 FloatT 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 
156 template<typename FloatT>
157 FloatT 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 
172 template<typename FloatT>
173 FloatT 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 
183 template<typename FloatT>
184 FloatT 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 
250 template<typename FloatT>
251 static 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 
268 template<typename FloatT>
269 static 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 
287 template<typename FloatT>
288 static 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 
305 template<typename FloatT>
306 static 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 
368 template
369 float GaussianConditional<float>(float mean, float stddev, float val);
370 template
371 float Gaussian<float>(float mean, float stddev, float val);
372 template
373 float PoissonConditional<float>(float mean, unsigned int step);
374 template
375 float Poisson<float>(float mean, unsigned int step);
376 template
377 float LogFactorial<float>(unsigned int n);
378 template
379 float Factorial<float>(unsigned int n);
380 template
381 float LogGamma<float>(float z);
382 template
383 float Gamma<float>(float z);
384 
385 template
386 double GaussianConditional<double>(double mean, double stddev, double val);
387 template
388 double Gaussian<double>(double mean, double stddev, double val);
389 template
390 double PoissonConditional<double>(double mean, unsigned int step);
391 template
392 double Poisson<double>(double mean, unsigned int step);
393 template
394 double LogFactorial<double>(unsigned int n);
395 template
396 double Factorial<double>(unsigned int n);
397 template
398 double LogGamma<double>(double z);
399 template
400 double 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.