218 lines
3.4 KiB
C++
Executable file
218 lines
3.4 KiB
C++
Executable file
//#include "StdAfx.h"
|
|
#include "Gamma.h"
|
|
#include "math.h"
|
|
#include <stdlib.h>
|
|
|
|
CGamma::CGamma()
|
|
{
|
|
|
|
}
|
|
|
|
CGamma::~CGamma()
|
|
{
|
|
|
|
}
|
|
|
|
double CGamma::gammln(double xx)
|
|
{
|
|
double RetVal;
|
|
|
|
double y;
|
|
double x;
|
|
double tmp;
|
|
double ser;
|
|
int j;
|
|
|
|
double cof[6];
|
|
|
|
cof[0] = 76.18009172947146; cof[1] = -86.50532032941677; cof[2] = 24.01409824083091;
|
|
cof[3] = -1.231739572450155; cof[4] = 1.208650973866179e-3; cof[5] = -5.395239384953e-6;
|
|
|
|
y = x = xx;
|
|
|
|
tmp = (x+5.5) - (x+0.5)*(log(x+5.5));
|
|
ser = 1.000000000190015;
|
|
|
|
for ( j=0; j<=5; j++ )
|
|
{
|
|
y += 1.0;
|
|
ser += cof[j]/y;
|
|
}
|
|
|
|
RetVal = log( 2.506628274631*ser/x ) - tmp;
|
|
|
|
return RetVal;
|
|
}
|
|
|
|
double CGamma::gamser(double a, double x)
|
|
{
|
|
double RetVal;
|
|
|
|
double gln;
|
|
double ap;
|
|
double del;
|
|
double sm;
|
|
int ITMAX;
|
|
int m;
|
|
|
|
if ( x!=0 )
|
|
{
|
|
ITMAX = 31;
|
|
gln = gammln( a );
|
|
ap = a;
|
|
del = sm = 1.0/a;
|
|
|
|
for( m=1; m<=ITMAX; m++ )
|
|
{
|
|
ap += 1;
|
|
del *= x/ap;
|
|
sm += del;
|
|
}
|
|
|
|
RetVal = sm * exp( a*log(x) - x - gln );
|
|
}
|
|
else
|
|
RetVal = 0;
|
|
|
|
return RetVal;
|
|
}
|
|
|
|
double CGamma::gcf(double a, double x)
|
|
{
|
|
double RetVal;
|
|
|
|
double gln;
|
|
double b;
|
|
double c;
|
|
double d;
|
|
double h;
|
|
double an;
|
|
double del;
|
|
double FPMIN;
|
|
int ITMAX;
|
|
int i;
|
|
|
|
gln = gammln( a );
|
|
b = x + 1.0 - a;
|
|
FPMIN = pow( 10., -30 );
|
|
c = 1.0 / FPMIN;
|
|
d = 1.0 / b;
|
|
h = d;
|
|
|
|
ITMAX = 30;
|
|
for ( i=1; i<=ITMAX; i++ )
|
|
{
|
|
an = (-i) * (i-a);
|
|
b += 2.0;
|
|
d = an*d + b;
|
|
if ( abs((int)d) < FPMIN ) d = FPMIN;
|
|
c = b + an/c;
|
|
if ( abs((int)c) < FPMIN ) c = FPMIN;
|
|
d = 1.0/d;
|
|
del = d * c;
|
|
h *= del;
|
|
}
|
|
|
|
RetVal = 1.0 - (exp( a*log(x) - x - gln ) * h);
|
|
|
|
return RetVal;
|
|
}
|
|
|
|
double CGamma::Gamma(double a, double chi2)
|
|
{
|
|
double RetVal;
|
|
|
|
double x;
|
|
|
|
a = a / 2;
|
|
x = chi2 / 2;
|
|
|
|
if ( x<=a )
|
|
RetVal = gamser( a, x );
|
|
else
|
|
RetVal = gcf( a, x );
|
|
|
|
return 1-RetVal;
|
|
}
|
|
|
|
double CGamma::Beta(double z, double a, double b)
|
|
{
|
|
double bt;
|
|
|
|
if (z < 0.0)
|
|
z = 0.0;
|
|
if (z > 1.0)
|
|
z = 1.0;
|
|
|
|
if (z == 0.0 || z == 1.0)
|
|
bt=0.0;
|
|
else
|
|
// Factors in front of the continued fraction.
|
|
bt = exp( gammln(a+b) - gammln(a) - gammln(b) + a*log(z) + b*log(1.0-z) );
|
|
|
|
if (z < (a+1.0)/(a+b+2.0)) // Use continued fraction directly.
|
|
return bt*betacf(a, b, z)/a;
|
|
else // Use continued fraction after making the symmetry transformation.
|
|
return 1.0-bt*betacf(b, a, 1.0-z)/b;
|
|
}
|
|
|
|
double CGamma::betacf(double a, double b, double x)
|
|
{
|
|
int m,m2;
|
|
double aa,c,d,del,h,qab,qam,qap;
|
|
|
|
double MAXIT = 100;
|
|
double EPS = 3.0e-7;
|
|
double FPMIN = 1.0e-30;
|
|
|
|
qab = a+b;
|
|
qap = a+1.0;
|
|
qam = a-1.0;
|
|
|
|
c = 1.0; // First step of Lentz's method.
|
|
d = 1.0 - qab*x/qap;
|
|
|
|
if (fabs(d) < FPMIN)
|
|
d = FPMIN;
|
|
|
|
d = 1.0/d;
|
|
h = d;
|
|
|
|
for (m=1; m<=MAXIT; m++)
|
|
// while (true)
|
|
{
|
|
m2 = 2*m;
|
|
aa = m*(b-m)*x / ((qam+m2)*(a+m2));
|
|
d = 1.0 + aa*d; // One step (the even one) of the recurrence.
|
|
|
|
if (fabs(d) < FPMIN)
|
|
d=FPMIN;
|
|
|
|
c = 1.0 + aa/c;
|
|
|
|
if (fabs(c) < FPMIN)
|
|
c = FPMIN;
|
|
|
|
d = 1.0/d;
|
|
h *= d*c;
|
|
aa = -(a+m)*(qab+m)*x / ((a+m2)*(qap+m2));
|
|
d = 1.0 + aa*d; // Next step of the recurrence (the odd one).
|
|
|
|
if (fabs(d) < FPMIN)
|
|
d = FPMIN;
|
|
|
|
c = 1.0 + aa/c;
|
|
|
|
if (fabs(c) < FPMIN)
|
|
c = FPMIN;
|
|
|
|
d = 1.0/d;
|
|
del = d*c;
|
|
h *= del;
|
|
|
|
if (fabs(del-1.0) < EPS) // Are we done?
|
|
break;
|
|
}
|
|
|
|
return h;
|
|
}
|