Файл gammaDF.h
/***********************************************************/
/* Gamma distribution */
/***********************************************************/
#ifndef __GAMMA_H__ /* To prevent redefinition */
#define ENTRY extern
#define LOCAL static
class GammaDF {
public:
GammaDF(double shape, double scale=1);
double value(double x); // Функция распределения Gamma(x|shape,scale)
double inv(double p); // Обратная функция: Gamma(x|shape,scale)=p
private:
double a, shape, scale, lga;
double fraction(double x);
double series(double x);
};
#define __GAMMA_H__ /* Prevents redefinition */
#endif /* Ends #ifndef __GAMMA_H__ */
Файл gammaDF.cpp
/***********************************************************/
/* Гамма-распределение */
/***********************************************************/
#include <math.h>
#include <assert.h>
#include "gammaDF.h"
#include "logGamma.h"
LOCAL const double zero = 0.0;
LOCAL const double one = 1.0;
LOCAL const double two = 2.0;
GammaDF::GammaDF(double shape, double scale):
a(shape), shape(shape), scale(scale), lga(logGamma(shape))
{
assert(shape > 0 && scale > 0);
}
double
GammaDF::series(double x)
// См. Abramowitz & Stegun,
// Handbook of Mathematical Functions, 1964 [6.5.29]
// М.Абрамовиц, И.Стиган
// Справочник по специальным функциям (М: Мир, 1979)
//
// Для вычисления неполной гамма функции P(a,x)
// используется ее разложение в ряд Тейлора.
//
{
double sum, prev_sum, term, aa = a;
long i = 0;
term = sum = one / a;
do {
aa += one;
term *= x / aa;
prev_sum = sum; sum += term;
i++;
} while(fabs(prev_sum) != fabs(sum));
return sum *= exp(-x + a * log(x) - lga);
}/* incGamma_series */
double
GammaDF::fraction(double x)
// См. Abramowitz & Stegun,
// Handbook of Mathematical Functions, 1964 [6.5.31]
// М.Абрамовиц, И.Стиган
// Справочник по специальным функциям (М: Мир, 1979)
//
// Для вычисления неполной гамма функции P(a,x)
// используется ее разложение в цепную дробь Лежандра
//
// P(a,x)=exp(-x +x*ln(a))*CF/logGamma(a),
//
// где
//
// 1 1-a 1 2-a 2
// CF = --- ----- --- ----- --- ....
// x+ 1+ x+ 1+ x+
//
// Используются подходящие дроби CF(n) = A(n) / B(n)
//
// где
// A(n) = (s(n) * A(n-1) + r(n) * A(n-2)) * factor
// B(n) = (s(n) * B(n-1) + r(n) * B(n-2)) * factor
// причем
// A(-1) = 1, B(-1) = 0, A(0) = s(0), B(0) = 1.
//
// Здесь
//
// s(0) = 0, s(1) = x, r(0) = 0, r(1) = 1,
//
// так что
//
// A(1) = one * factor, B(1) = r * factor
//
// и, следовательно,
//
// r(i) = k - a if i = 2k, k > 0
// r(i) = k if i = 2k+1,
// s(i) = 1 if i = 2k,
// s(i) = x if i = 2k+1
//
// factor - шкалирующий множитель
//
{
double old_sum=zero, factor=one;
double A0=zero, A1=one, B0=one, B1=x;
double sum=one/x, z=zero, ma=zero-a, rfact;
do {
z += one;
ma += one;
/* two steps of recurrence replace A's & B's */
A0 = (A1 + ma * A0) * factor; /* i even */
B0 = (B1 + ma * B0) * factor;
rfact = z * factor;
A1 = x * A0 + rfact * A1; /* i odd, A0 already rescaled */
B1 = x * B0 + rfact * B1;
if (B1) {
factor = one / B1;
old_sum = sum;
sum = A1 * factor;
}/* if */
} while (fabs(sum) != fabs(old_sum));
return exp(-x + a * log(x) - lga) * sum;
}/*fraction*/
double
GammaDF::value(double x)
// Вычисляется Gamma(x|a):
// вероятность того, что случайная величина,
// подчиняющаяся центральному гамма-распределению с параметром 'a',
// меньше или равна 'x'.
{
x *= scale;
if(x <= zero)
return zero; /* НЕ ошибка! */
else if(x < (a + one))
return series(x);
else
return one - fraction(x);
}/*value*/
double
GammaDF::inv(double p)
// Ищет такое значение 'x', для которого Gamma(x|a) = p,
// т.е. равна 'p' вероятность того, что случайная величина,
// подчиняющаяся центральному гамма-распределению с параметром 'a',
// меньше или равна 'x'.
{
double fx, l = 0, r = 1, x;
if (p == 0) return 0;
assert(p > 0 && p < 1);
for(l=0, r=a/2; value(r) < p; r+=a/2) l=r;
x=(l+r)/2;
do {
fx = value(x);
if (fx > p) r = x; else
if (fx < p) l = x; else
break;
x = (l + r)* 0.5;
} while ((l!=x) && (r!=x));
return x;
}/*inv*/
1 2
8 8 8
| 
*
Обратная связь