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);
Generic library namespace.
FloatT Gamma(FloatT z)
Euler's Gamma function.
static FloatType sqrt_pi()
The square root of pi.
FloatT Factorial(unsigned int n)
Gives n!
static FloatType sqrt2()
The square root of 2.
FloatT GaussianConditional(FloatT mean, FloatT stddev, FloatT val)
Gives the conditional probability of the Gaussian distribution at position val.
static FloatType pi()
The constant pi.
FloatT LogFactorial(unsigned int n)
Gives the natural log of n!
FloatT Poisson(FloatT mean, unsigned int step)
Gives the value of the Poisson distribution at position step.
static FloatType log2()
The natural logarithm of 2.
static FloatType log_pi()
The natural logarithm of pi.
FloatT Gaussian(FloatT mean, FloatT stddev, FloatT val)
Gives the value of the Gaussian distribution at position val.
FloatT PoissonConditional(FloatT mean, unsigned int step)
Gives the conditional probability of the Poisson distribution at position step.
FloatT LogGamma(FloatT z)
The natural log of Euler's Gamma function.