26 #include <wfmath/probability.h>
28 #include <wfmath/const.h>
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);
48 static const unsigned int GammaCutoff = 10;
50 template<
typename FloatT>
55 FloatT diff = val - mean;
56 FloatT diffnorm = diff / stddev;
57 FloatT diffsqr_over_two = diffnorm * diffnorm / 2;
62 static const FloatT half = 0.5;
64 10 * std::numeric_limits<FloatT>::epsilon()) {
65 FloatT erfc_norm = IncompleteGammaComplement(half, diffsqr_over_two);
67 FloatT normalization = (diffnorm > 0) ? (erfc_norm / 2) : (1 - erfc_norm / 2);
69 return Gaussian(mean, stddev, val) / normalization;
71 static const FloatT two = 2.0;
73 return two / (std::fabs(diff)
74 * IncompleteGammaComplementNoPrefactor<FloatT>(half, diffsqr_over_two));
77 template<
typename FloatT>
78 FloatT
Gaussian(FloatT mean, FloatT stddev, FloatT val)
82 FloatT diff = (mean - val) / stddev;
84 return std::exp(-(diff * diff) / 2) / (std::fabs(stddev) *
89 template<
typename FloatT>
95 return (step == 0) ? 1 : 0;
98 return std::exp(-mean);
102 IncompleteGamma(
static_cast<FloatT
>(step), mean);
104 static const FloatT one = 1.0;
105 return one / IncompleteGammaNoPrefactor(
static_cast<FloatT
>(step), mean);
108 template<
typename FloatT>
109 FloatT
Poisson(FloatT mean,
unsigned int step)
114 return (step == 0) ? 1 : 0;
116 return std::exp(LogPoisson(mean, step));
119 template<
typename FloatT>
120 static FloatT LogPoisson(FloatT mean,
unsigned int step)
127 FloatT first =
static_cast<FloatT
>(step) * std::log(mean);
128 FloatT second = mean + LogFactorial<FloatT>(step);
130 assert(
"LogFactorial() always returns positive" && second > 0);
132 FloatT ans = first - second;
135 assert(
"must get probability < 1" && ans < 0);
140 template<
typename FloatT>
146 if(n < GammaCutoff) {
147 FloatT ans =
static_cast<FloatT
>(n);
149 ans *=
static_cast<FloatT
>(n);
153 return std::exp(
LogGamma(
static_cast<FloatT
>(n + 1)));
156 template<
typename FloatT>
162 if(n < GammaCutoff) {
163 FloatT ans =
static_cast<FloatT
>(n);
165 ans *=
static_cast<FloatT
>(n);
166 return std::log(ans);
169 return LogGamma(
static_cast<FloatT
>(n + 1));
172 template<
typename FloatT>
183 template<
typename FloatT>
190 static const FloatT one = 1.0;
194 LogGamma<FloatT>(one - z) -
205 if(z < GammaCutoff) {
207 while(z < GammaCutoff)
209 log_shift = std::log(std::fabs(shift));
217 static const FloatT half = 0.5, two = 2.0;
219 FloatT ans = (z - half) * std::log(z) - log_shift - z +
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);
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();
235 for(i = 0; i < num_coeffs; ++i) {
236 FloatT next_term = coeffs[i] * z_power;
238 if(std::fabs(next_term) < small_enough)
240 z_power *= z_to_minus_two;
243 assert(i < num_coeffs);
250 template<
typename FloatT>
251 static FloatT IncompleteGamma(FloatT a, FloatT z)
253 assert(z >= 0 && a >= 0);
261 return 1 - IncompleteGammaComplement(a, z);
263 FloatT prefactor = std::exp(a * (std::log(z) + 1) - z -
LogGamma(a));
265 return IncompleteGammaNoPrefactor(a, z) * prefactor;
268 template<
typename FloatT>
269 static FloatT IncompleteGammaNoPrefactor(FloatT a, FloatT z)
271 assert(z <= a + 1 && z > 0 && a > 0);
279 while(std::fabs(term / ans) > std::numeric_limits<FloatT>::epsilon()) {
280 term *= z / ++dividend;
287 template<
typename FloatT>
288 static FloatT IncompleteGammaComplement(FloatT a, FloatT z)
290 assert(z >= 0 && a >= 0);
298 return 1 - IncompleteGamma(a, z);
300 FloatT prefactor = std::exp(a * std::log(z) - z -
LogGamma(a));
302 return IncompleteGammaComplementNoPrefactor(a, z) * prefactor;
305 template<
typename FloatT>
306 static FloatT IncompleteGammaComplementNoPrefactor(FloatT a, FloatT z)
308 assert(z >= a + 1 && a > 0);
312 static const FloatT fudge = 1000;
314 FloatT b_contrib = z + 1 - a;
315 FloatT a_last, b_last, a_next, b_next;
318 next_zero = (std::fabs(b_contrib) <= std::numeric_limits<FloatT>::min()
329 b_last = 1 / b_contrib;
335 FloatT a_tmp = a_next, b_tmp = b_next;
337 FloatT a_contrib = term * (a - term);
341 a_next = b_contrib * a_tmp + a_contrib * a_last;
342 b_next = b_contrib * b_tmp + a_contrib * b_last;
347 last_zero = next_zero;
348 next_zero = (std::fabs(b_next) <=
349 std::fabs(a_next) * (std::numeric_limits<FloatT>::min() *
358 std::fabs(a_next - a_last) < std::fabs(a_last) *
359 std::numeric_limits<FloatT>::epsilon())
369 float GaussianConditional<float>(
float mean,
float stddev,
float val);
371 float Gaussian<float>(
float mean,
float stddev,
float val);
373 float PoissonConditional<float>(
float mean,
unsigned int step);
375 float Poisson<float>(
float mean,
unsigned int step);
377 float LogFactorial<float>(
unsigned int n);
379 float Factorial<float>(
unsigned int n);
381 float LogGamma<float>(
float z);
383 float Gamma<float>(
float z);
386 double GaussianConditional<double>(
double mean,
double stddev,
double val);
388 double Gaussian<double>(
double mean,
double stddev,
double val);
390 double PoissonConditional<double>(
double mean,
unsigned int step);
392 double Poisson<double>(
double mean,
unsigned int step);
394 double LogFactorial<double>(
unsigned int n);
396 double Factorial<double>(
unsigned int n);
398 double LogGamma<double>(
double z);
400 double Gamma<double>(
double z);