    001    package calhoun.analysis.crf.statistics;
    002
    003    import calhoun.util.Assert;
    004
    005    public class GammaDistribution {
    006
    007            /* See http://mathworld.wolfram.com/GammaDistribution.html
    008             *
    009             *
    010             *
    011             *
    012             */
    013
    014            public static double lgamma(double p, double lambda, double x) {
    015                    Assert.a(lambda>0);
    016                    Assert.a(x>0);
    017                    Assert.a(p>0);
    018                    double ret = (p-1)*Math.log(x) -x*lambda + p*Math.log(lambda) - lgamma(p);
    019                    if (Double.isNaN(ret)) {
    020                            Assert.a(false,"  pow1=" +Math.pow(x,p-1)+"  exp1="+Math.exp(-x*lambda)+"  pow2="+Math.pow(lambda,p) + "  gamma="+gamma(p) );
    021                    }
    022                    return ret;
    023            }
    024
    025
    026            public static double gamma(double p, double lambda, double x) {
    027                    Assert.a(lambda>0);
    028                    Assert.a(x>0);
    029                    Assert.a(p>0);
    030
    031                    double lret = lgamma(p,lambda,x);
    032                    double ret = Math.exp(lret);
    033                    //Assert.a(!Double.isNaN(ret));
    034                    Assert.a((ret != Double.NEGATIVE_INFINITY) && (ret != Double.POSITIVE_INFINITY) && (!Double.isNaN(ret)));
    035                    return ret;
    036            }
    037
    038            public static double gamma(double x) {
    039                    /* The gamma function, which is the improper integral
    040                     * integral_0^infinity t^(x-1) * exp(-t) dt
    041                     *
    042                     *      x*gamma(x) = gamma(x+1).
    043                     *
    044                     * For integers x>=1, gamma(x) = (x-1)!
    045                     */
    046                    Assert.a(x>0);
    047                    double y = lgamma(x);
    048                    return Math.exp(y);
    049            }
    050
    051
    052             private static double lgamma(double x) {
    053                     Assert.a(x>0);
    054
    055                     double l2pi = Math.log(2*Math.PI);
    056                     double term, y, zz ;
    057                     int k, t ;
    058
    059                     if (x<10.0) return (lgamma(x+1.0) - Math.log(x)) ;
    060                     y = (x-0.5)*Math.log(x) -x + 0.5*l2pi ;
    061                     zz = 1.0/x ;
    062                     for (k=1; (2*k)<= bern.length ; k++)  {
    063                      t = 2*k ;
    064                      term = bern[2*k]/(double) (t*(t-1)) ;
    065                      term *= zz ;
    066                      y += term ;
    067                      zz /= (x*x) ;
    068                     }
    069
    070                    return y;
    071            }
    072
    073
    074
    075
    076            // Ported from Nick Patterson's code 20060203
    077             static private double[] bern = new double[]{1.0,
    078                             -1.0/2.0,
    079                             1.0/6.0,
    080                             0,
    081                             -1.0/30.0,
    082                             0,
    083                             1.0/42.0,
    084                             0,
    085                             -1.0/30.0,
    086                             0,
    087                             5.0/66.0,
    088                             0,
    089                             -691.0/2730.0,
    090                             0,
    091                             7.0/6.0
    092             };
    093
    094             // Ported from Nick Patterson's code 20060203
    095             public static double[] mleg(double a1, double a2) // both p and lam assumed of size 1
    096            //  solve
    097            //  p/lam = a1 ; psi(p) - log(lam) = a2 ;
    098            //  Thus psi(p) - log(p) = a2 - log(a1)
    099            {
    100               int iter ;
    101               double s, pp, ll  ;
    102               double top, bot, fval ;
    103
    104
    105              s = a2 - Math.log(a1) ;
    106
    107              Assert.a(s<=0.0 , "log E(x) < E(log (x)) \n" ) ;
    108              pp = -s ;
    109
    110              for (iter = 1; iter <= 30 ; ++iter) {
    111               fval = s - (psi(pp) - Math.log (pp)) ;
    112               if (fval<0.0)  break ;
    113               pp *= 2.0 ;
    114              }
    115
    116              for (iter = 1; iter <= 30 ; ++iter) {
    117               fval = s - (psi(pp) - Math.log (pp)) ;
    118               if (fval>0.0)  break ;
    119               pp /= 2.0 ;
    120              }
    121
    122              for (iter = 1; iter <= 10 ; ++iter) {
    123               fval = psi(pp) - Math.log (pp) ;
    124               top = s-fval ;
    125               bot =  tau(pp) - (1.0/pp) ;
    126               pp += top/bot ;
    127              }
    128              ll = pp/a1 ;
    129
    130              return (new double[]{pp,ll});
    131            }
    132
    133
    134
    135            private static double psi(double x) {
    136                    // Ported from Nick Pattersons code 20060203
    137                    // psi is the derivative of the gamma function
    138
    139                     double y, zz, term ;
    140                     int k ;
    141
    142                     Assert.a( x > 0.0 , "(psi) bad value: " +  x) ;
    143
    144                     Assert.a(bern.length == 15);
    145
    146                     if (x<10.0) return (psi(x+1.0) - 1.0/x) ;
    147
    148                     y = Math.log(x) - 1.0/(2.0*x) ;
    149                     zz = 1.0 ;
    150                     for (k=1; (2*k)< bern.length ; k++)  {
    151                      zz /= (x*x) ;
    152                      term = bern[2*k]/(double) (2*k) ;
    153                      term *= zz ;
    154                      y -= term ;
    155                     }
    156                     return y ;
    157            }
    158
    159
    160            private static double tau(double x)
    161            //       Ported from Nick Pattersons code 20060203
    162            //      derivative of psi
    163            {
    164                     double y, zz, term ;
    165                     int k ;
    166
    167                     Assert.a(x>0.0 , "(tau) bad value: " +  x) ;
    168
    169                     if (x<10.0) return (tau(x+1.0) + 1.0/(x*x)) ;
    170
    171                     y = 1.0/x  + 1.0/(2.0*x*x) ;
    172                     zz = 1.0/x ;
    173                     for (k=1; (2*k)< bern.length ; k++)  {
    174                      zz /= (x*x) ;
    175                      term = bern[2*k]/(double) (2*k) ;
    176                      term *= zz ;
    177                      term *= - (double) (2*k) ;
    178                      y -= term ;
    179                     }
    180                     return y ;
    181            }
    182
    183
    184    }