mirror of
https://github.com/UpsilonNumworks/Upsilon.git
synced 2026-01-19 08:47:28 +01:00
82 lines
3.2 KiB
C++
82 lines
3.2 KiB
C++
#include "regularized_gamma.h"
|
|
#include "distribution.h"
|
|
#include "helper.h"
|
|
#include <cmath>
|
|
#include <float.h>
|
|
#include <assert.h>
|
|
|
|
bool regularizedGamma(double s, double x, double epsilon, int maxNumberOfIterations, double * result) {
|
|
// TODO Put interruption instead of maxNumberOfIterations
|
|
|
|
assert(!std::isnan(s) && !std::isnan(x) && s > 0.0 && x >= 0.0);
|
|
if (x == 0.0) {
|
|
*result = 0.0;
|
|
return true;
|
|
}
|
|
if (std::isinf(x)) {
|
|
*result = 1.0;
|
|
return true;
|
|
}
|
|
if (x >= s + 1.0) {
|
|
/* For x >= s + 1: Compute using the continued fraction representation
|
|
* regularizedGamma(s,x) = 1 - gamma(s,x)/gamma(s)
|
|
* 1
|
|
* = 1 - e^-x * x^s / gamma(s) * -----------------------
|
|
* 1(1-s)
|
|
* 1+x-s- ----------------
|
|
* 2(2-s)
|
|
* 3+x-s- ---------
|
|
* x+5-s-...
|
|
*
|
|
* 1
|
|
* = 1 - e^-x * x^s / gamma(s) * -----------------------
|
|
* a(1)
|
|
* b(0)- ----------------
|
|
* a(2)
|
|
* b(1)- ---------
|
|
* b(3)-...
|
|
*
|
|
* With a(n) = n(n-s)
|
|
* With b(n) = (2*n+1) + x - s
|
|
*
|
|
* We have the following formulae:
|
|
* A(n)=b(n)A(n-1)+a(n)A(n-2) A(-1) = 1 A(0) = b(0)
|
|
* B(n)=b(n)B(n-1)+a(n)B(n-2) B(-1) = 0 B(0) = 1
|
|
* The consecutive convergents of the continued fraction are : A(n)/B(n)
|
|
* => This method does not work well.
|
|
*
|
|
* Other method : Lentz's method with Thompson and Barnett's modification */
|
|
|
|
double continuedFractionValue = 0.0;
|
|
if (!Helper::ContinuedFractionEvaluation(
|
|
[](double index, double s, double x) { return index * (s - index); },
|
|
[](double index, double s, double x) { return (2.0 * index + 1.0) + x - s; },
|
|
maxNumberOfIterations,
|
|
&continuedFractionValue,
|
|
s,
|
|
x))
|
|
{
|
|
return false;
|
|
}
|
|
*result = 1.0 - std::exp(-x + s*std::log(x) - std::lgamma(s)) * ( 1.0 / continuedFractionValue);
|
|
return true;
|
|
}
|
|
|
|
/* For x < s + 1: Compute using the infinite series representation
|
|
* incompleteGamma(s,x) = 1/gamma(a) * e^-x * x^s * sum(x^n/((s+n)(s+n-1)...) */
|
|
double infiniteSeriesValue = 0.0;
|
|
if (!Helper::InfiniteSeriesEvaluation(
|
|
1.0/s,
|
|
[](double previousTerm, double index, double s, double x, double d1, double d2) { return previousTerm * x / (s + index); },
|
|
epsilon,
|
|
maxNumberOfIterations,
|
|
&infiniteSeriesValue,
|
|
s,
|
|
x))
|
|
{
|
|
return false;
|
|
}
|
|
*result = std::isinf(infiniteSeriesValue) ? 1.0 : std::exp(-x + s*std::log(x) - std::lgamma(s)) * infiniteSeriesValue;
|
|
return true;
|
|
}
|