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
WFMath::LogFactorial
FloatT LogFactorial(unsigned int n)
Gives the natural log of n!
Definition: probability.cpp:157
WFMath::Poisson
FloatT Poisson(FloatT mean, unsigned int step)
Gives the value of the Poisson distribution at position step.
Definition: probability.cpp:109
WFMath::numeric_constants::sqrt_pi
static FloatType sqrt_pi()
The square root of pi.
WFMath::numeric_constants::sqrt2
static FloatType sqrt2()
The square root of 2.
WFMath::Gaussian
FloatT Gaussian(FloatT mean, FloatT stddev, FloatT val)
Gives the value of the Gaussian distribution at position val.
Definition: probability.cpp:78
WFMath::PoissonConditional
FloatT PoissonConditional(FloatT mean, unsigned int step)
Gives the conditional probability of the Poisson distribution at position step.
Definition: probability.cpp:90
WFMath::numeric_constants::pi
static FloatType pi()
The constant pi.
WFMath::LogGamma
FloatT LogGamma(FloatT z)
The natural log of Euler's Gamma function.
Definition: probability.cpp:184
WFMath
Generic library namespace.
Definition: shape.h:41
WFMath::numeric_constants::log2
static FloatType log2()
The natural logarithm of 2.
WFMath::Gamma
FloatT Gamma(FloatT z)
Euler's Gamma function.
Definition: probability.cpp:173
WFMath::numeric_constants
Definition: const.h:65
WFMath::Factorial
FloatT Factorial(unsigned int n)
Gives n!
Definition: probability.cpp:141
WFMath::numeric_constants::log_pi
static FloatType log_pi()
The natural logarithm of pi.
WFMath::GaussianConditional
FloatT GaussianConditional(FloatT mean, FloatT stddev, FloatT val)
Gives the conditional probability of the Gaussian distribution at position val.
Definition: probability.cpp:51