root/trunk/mathextra/gammadist.d

// incomplete gamma functions

module mathextra.gammadist;
private import std.math;
private import mathextra.normaldist;
debug import std.stdio;

unittest {
    debug writefln("--- UnitTest: " __FILE__ " ---");
} 

// There are still problems with std.math.tgamma.
private import mathextra.etcgamma;
alias mathextra.etcgamma.tgamma gamma;
alias mathextra.etcgamma.lgamma loggamma;


real gammaIncomplete(real a, real x )
in {
   assert(x >= 0);
   assert(a > 0);
}
body {
    /* left tail of incomplete gamma function:
     *
     *          inf.      k
     *   a  -x   -       x
     *  x  e     >   ----------
     *           -     -
     *          k=0   | (a+k+1)
     *
     */
    if (x==0)
       return 0.0L;

    if ( (x > 1.0L) && (x > a ) )
        return 1.0L - gammaIncompleteCompl(a,x);
    
    real ax = a * log(x) - x - loggamma(a);
/+  
    if( ax < MINLOGL ) return 0; // underflow
    //  { mtherr( "igaml", UNDERFLOW ); return( 0.0L ); }
+/  
    ax = exp(ax);
    
    /* power series */
    real r = a;
    real c = 1.0L;
    real ans = 1.0L;
    
    do  {
        r += 1.0L;
        c *= x/r;
        ans += c;
    } while( c/ans > real.epsilon );
    
    return ans * ax/a;
}

/** ditto */
real gammaIncompleteCompl(real a, real x )
in {
   assert(x >= 0);
   assert(a > 0);
}
body {
    if (x==0)
       return 1.0L;
    if ( (x < 1.0L) || (x < a) )
        return 1.0L - gammaIncomplete(a,x);

   // DAC (Cephes bug fix): This is necessary to avoid
   // spurious nans, eg
   // log(x)-x = NaN when x = real.infinity
    const real MAXLOGL =  1.1356523406294143949492E4L;
   if (x > MAXLOGL) return 0; // underflow
    
    real ax = a * log(x) - x - loggamma(a);
//const real MINLOGL = -1.1355137111933024058873E4L;
//  if ( ax < MINLOGL ) return 0; // underflow;
    ax = exp(ax);

    
    /* continued fraction */
    real y = 1.0L - a;
    real z = x + y + 1.0L;
    real c = 0.0L;
    
    real pk, qk, t;

    real pkm2 = 1.0L;
    real qkm2 = x;
    real pkm1 = x + 1.0L;
    real qkm1 = z * x;
    real ans = pkm1/qkm1;
    
    do  {   
        c += 1.0L;
        y += 1.0L;
        z += 2.0L;
        real yc = y * c;
        pk = pkm1 * z  -  pkm2 * yc;
        qk = qkm1 * z  -  qkm2 * yc;
        if( qk != 0.0L ) {
            real r = pk/qk;
            t = fabs( (ans - r)/r );
            ans = r;
        } else {
            t = 1.0L;
        }
    pkm2 = pkm1;
        pkm1 = pk;
        qkm2 = qkm1;
        qkm1 = qk;
        
        const real BIG = 9.223372036854775808e18L;
        
        if ( fabs(pk) > BIG ) {
            pkm2 /= BIG;
            pkm1 /= BIG;
            qkm2 /= BIG;
            qkm1 /= BIG;
        }
    } while ( t > real.epsilon );
    
    return ans * ax;
}

/*
 *      Inverse of complemented incomplete gamma integral
 *
 * Given a and y, the function finds x such that
 *
 *  gammaIncompleteCompl( a, x ) = p.
 *
 * Starting with the approximate value x = a $(POWER t, 3), where
 * t = 1 - d - normalDistributionInv(p) sqrt(d), 
 * and d = 1/9a,
 * the routine performs up to 10 Newton iterations to find the
 * root of incompleteGammaCompl(a,x) - p = 0.
 */

// MINLOGL = -1.1355137111933024058873E4L;
real gammaIncompleteComplInv(real a, real p)
in {
  assert(p>=0 && p<= 1);
  assert(a>0);
}
body {
    if (p==0) return real.infinity;
    
    real y0 = p;
    const real MAXLOGL =  1.1356523406294143949492E4L;
    real x0, x1, x, yl, yh, y, d, lgm, dithresh;
    int i, dir;
    
    /* bound the solution */
    x0 = real.max;
    yl = 0.0L;
    x1 = 0.0L;
    yh = 1.0L;
    dithresh = 4.0 * real.epsilon;
    
    /* approximation to inverse function */
    d = 1.0L/(9.0L*a);
    y = 1.0L - d - normalDistributionInv(y0) * sqrt(d);
    x = a * y * y * y;
    
    lgm = loggamma(a);
    
    for( i=0; i<10; i++ ) {
        if( x > x0 || x < x1 )
            goto ihalve;
        y = gammaIncompleteCompl(a,x);
        if ( y < yl || y > yh )
            goto ihalve;
        if ( y < y0 ) {
            x0 = x;
            yl = y;
        } else {
            x1 = x;
            yh = y;
        }
    /* compute the derivative of the function at this point */
        d = (a - 1.0L) * log(x0) - x0 - lgm;
        if ( d < -MAXLOGL )
            goto ihalve;
        d = -exp(d);
    /* compute the step to the next approximation of x */
        d = (y - y0)/d;
        x = x - d;
        if ( i < 3 ) continue;
        if ( fabs(d/x) < dithresh ) return x;
    }
    
    /* Resort to interval halving if Newton iteration did not converge. */
ihalve:
    d = 0.0625L;
    if ( x0 == real.max ) {
        if( x <= 0.0L )
            x = 1.0L;
        while( x0 == real.max ) {
            x = (1.0L + d) * x;
            y = gammaIncompleteCompl( a, x );
            if ( y < y0 ) {
                x0 = x;
                yl = y;
                break;
            }
            d = d + d;
        }
    }
    d = 0.5L;
    dir = 0;

    for( i=0; i<400; i++ ) {
        x = x1  +  d * (x0 - x1);
        y = gammaIncompleteCompl( a, x );
        lgm = (x0 - x1)/(x1 + x0);
        if ( fabs(lgm) < dithresh )
            break;
        lgm = (y - y0)/y0;
        if ( fabs(lgm) < dithresh )
            break;
        if ( x <= 0.0L )
            break;
        if ( y > y0 ) {
            x1 = x;
            yh = y;
            if ( dir < 0 ) {
                dir = 0;
                d = 0.5L;
            } else if ( dir > 1 )
                d = 0.5L * d + 0.5L; 
            else
                d = (y0 - yl)/(yh - yl);
            dir += 1;
        } else {
            x0 = x;
            yl = y;
            if ( dir > 0 ) {
                dir = 0;
                d = 0.5L;
            } else if ( dir < -1 )
                d = 0.5L * d;
            else
                d = (y0 - yl)/(yh - yl);
            dir -= 1;
        }
    }
    /+
    if( x == 0.0L )
        mtherr( "igamil", UNDERFLOW );
    +/
    return x;
}

unittest {
//Values from Excel's GammaInv(1-p, x, 1)
assert(fabs(gammaIncompleteComplInv(1, 0.5) - 0.693147188044814) < 0.00000005);
assert(fabs(gammaIncompleteComplInv(12, 0.99) - 5.42818075054289) < 0.00000005);
assert(fabs(gammaIncompleteComplInv(100, 0.8) - 91.5013985848288L) < 0.000005);

assert(gammaIncomplete(1, 0)==0);
assert(gammaIncompleteCompl(1, 0)==1);
assert(gammaIncomplete(4545, real.infinity)==1);

// Values from Excel's (1-GammaDist(x, alpha, 1, TRUE))

assert(fabs(1.0L-gammaIncompleteCompl(0.5, 2) - 0.954499729507309L) < 0.00000005);
assert(fabs(gammaIncomplete(0.5, 2) - 0.954499729507309L) < 0.00000005);
// Fixed Cephes bug:
assert(gammaIncompleteCompl(384, real.infinity)==0);
assert(gammaIncompleteComplInv(3, 0)==real.infinity);

//writefln("%.20g",gammaIncompleteCompl(100, 0));
//assert(gammaIncompleteComplInv(8, 0));

// BUG: infinite loop if p == 0!
//writefln(gammaIncompleteComplInv(8, 0));

//writefln(gammaIncompleteComplInv(8, 1e-50));
//writefln(gammaIncompleteComplInv(12, 0.99));

}