26#include <wfmath/probability.h>
28#include <wfmath/const.h>
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);
48static const unsigned int GammaCutoff = 10;
50template<
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));
77template<
typename FloatT>
78FloatT
Gaussian(FloatT mean, FloatT stddev, FloatT val)
82 FloatT diff = (mean - val) / stddev;
84 return std::exp(-(diff * diff) / 2) / (std::fabs(stddev) *
89template<
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);
108template<
typename FloatT>
114 return (step == 0) ? 1 : 0;
116 return std::exp(LogPoisson(mean, step));
119template<
typename FloatT>
120static 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);
140template<
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)));
156template<
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));
172template<
typename FloatT>
183template<
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);
250template<
typename FloatT>
251static 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;
268template<
typename FloatT>
269static 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;
287template<
typename FloatT>
288static 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;
305template<
typename FloatT>
306static 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())
369float GaussianConditional<float>(
float mean,
float stddev,
float val);
371float Gaussian<float>(
float mean,
float stddev,
float val);
373float PoissonConditional<float>(
float mean,
unsigned int step);
375float Poisson<float>(
float mean,
unsigned int step);
377float LogFactorial<float>(
unsigned int n);
379float Factorial<float>(
unsigned int n);
381float LogGamma<float>(
float z);
383float Gamma<float>(
float z);
386double GaussianConditional<double>(
double mean,
double stddev,
double val);
388double Gaussian<double>(
double mean,
double stddev,
double val);
390double PoissonConditional<double>(
double mean,
unsigned int step);
392double Poisson<double>(
double mean,
unsigned int step);
394double LogFactorial<double>(
unsigned int n);
396double Factorial<double>(
unsigned int n);
398double LogGamma<double>(
double z);
400double Gamma<double>(
double z);
Generic library namespace.
FloatT GaussianConditional(FloatT mean, FloatT stddev, FloatT val)
Gives the conditional probability of the Gaussian distribution at position val.
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.
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.
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.