/***********************************************************/
/*                   Гамма-распределение                   */
/***********************************************************/

#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*/