ADMB Documentation  11.1.2246
 All Classes Files Functions Variables Typedefs Friends Defines
cgamdev.cpp
Go to the documentation of this file.
00001 
00016 #include <fvar.hpp>
00017 #define ITMAX 100
00018 #define EPS 1.e-9
00019 //#define EPS 3.0e-7
00020 #define FPMIN 1.0e-30
00021 void gcf(double& gammcf,double a,double x,double &gln);
00022 void gser(double& gamser,double a,double x,double& gln);
00023 
00024 double gammp(double a,double x)
00025 {
00026   double gamser,gammcf,gln;
00027 
00028   if (x < 0.0 || a <= 0.0)
00029     cerr << "Invalid arguments in routine gammp" << endl;
00030   if (x < (a+1.0)) {
00031     gser(gamser,a,x,gln);
00032     return gamser;
00033   } else {
00034     gcf(gammcf,a,x,gln);
00035     return 1.0-gammcf;
00036   }
00037 }
00038 
00039 double cumd_gamma(double x,double a)
00040 {
00041   double gamser,gammcf,gln;
00042 
00043   if (x < 0.0 || a <= 0.0)
00044     cerr << "Invalid arguments in routine gammp" << endl;
00045   if (x < (a+1.0)) {
00046     gser(gamser,a,x,gln);
00047     return gamser;
00048   } else {
00049     gcf(gammcf,a,x,gln);
00050     return 1.0-gammcf;
00051   }
00052 }
00053 
00060 void gcf(double& gammcf,double a,double x,double &gln)
00061 {
00062   int i;
00063   double an,b,c,d,del,h;
00064 
00065   gln=gammln(a);
00066   b=x+1.0-a;
00067   c=1.0/FPMIN;
00068   d=1.0/b;
00069   h=d;
00070   for (i=1;i<=ITMAX;i++) {
00071     an = -i*(i-a);
00072     b += 2.0;
00073     d=an*d+b;
00074     if (fabs(d) < FPMIN) d=FPMIN;
00075     c=b+an/c;
00076     if (fabs(c) < FPMIN) c=FPMIN;
00077     d=1.0/d;
00078     del=d*c;
00079     h *= del;
00080     if (fabs(del-1.0) < EPS) break;
00081   }
00082   if (i > ITMAX)
00083     cerr << "a too large, ITMAX too small in gcf" << endl;
00084   gammcf=exp(-x+a*log(x)-(gln))*h;
00085 }
00086 
00093 void gser(double& gamser,double a,double x,double& gln)
00094 {
00095   int n;
00096   double sum,del,ap;
00097 
00098   gln=gammln(a);
00099   if (x <= 0.0) {
00100     if (x < 0.0)
00101       cerr << "x less than 0 in routine gser" << endl;
00102     gamser=0.0;
00103     return;
00104   } else {
00105     ap=a;
00106     del=sum=1.0/a;
00107     for (n=1;n<=ITMAX;n++) {
00108       ++ap;
00109       del *= x/ap;
00110       sum += del;
00111       if (fabs(del) < fabs(sum)*EPS) {
00112         gamser=sum*exp(-x+a*log(x)-(gln));
00113         return;
00114       }
00115     }
00116     cerr << "a too large, ITMAX too small in routine gser" << endl;
00117     return;
00118   }
00119 }
00120 
00121 double get_initial_u(double a,double y);
00122 
00123 double Sn(double x,double a);
00124 
00125 double inv_cumd_gamma(double y,double a)
00126 {
00127   if (a<0.05)
00128   {
00129     cerr << "a musdt be > 0.1" << endl;
00130     ad_exit(1);
00131   }
00132   double u=get_initial_u(a,y);
00133   double h;
00134   do
00135   {
00136     double z=gammp(a,a*exp(u));
00137     double d=y-z;
00138     //cout << d << endl;
00139     double log_fprime=a*log(a)+a*(u-exp(u)) -gammln(a);
00140     double fprime=exp(log_fprime);
00141     h=d/fprime;
00142     u+=h;
00143   }
00144   while(fabs(h)>1.e-12);
00145   double x=a*exp(u);
00146   return x;
00147 }
00148 
00149 #undef ITMAX
00150 #undef EPS
00151 
00152 double Sn(double x,double a)
00153 {
00154   int i=1;
00155   double xp=x;
00156   double prod=1.0;
00157   double summ=1.0;
00158   double summand;
00159   do
00160   {
00161     prod*=(a+i);
00162     summand=xp/prod;
00163     if (summand<1.e-4) break;
00164     summ+=summand;
00165     i++;
00166     if (i>50)
00167     {
00168       cerr << "convergence error" << endl;
00169       ad_exit(1);
00170     }
00171   }
00172   while(1);
00173   return summ;
00174 }
00175 
00176 double get_initial_u(double a,double y)
00177 {
00178   const double c=0.57721;
00179   // note that P = y;
00180   double logP=log(y);
00181   double logQ=log(1-y);
00182   double logB=logQ+gammln(a);
00183   double x0=1.e+100;
00184   double log_x0=1.e+100;
00185 
00186   if (a<1.0)
00187   {
00188     if ( logB>log(.6) || (logB > log(.45) && a>=.3) )
00189     {
00190       double logu;
00191       if (logB+logQ > log(1.e-8))
00192       {
00193         logu=(logP+gammln(1.0+a))/a;
00194       }
00195       else
00196       {
00197         logu=-exp(logQ)/a -c;
00198       }
00199       double u=exp(logu);
00200       x0=u/(1-u/(1.0+a));
00201       double tmp=log(1-u/(1.0+a));
00202       log_x0=logu;
00203       log_x0-=tmp;
00204     }
00205     else if ( a<.3 && log(.35) <= logB && logB <= log(.6) )
00206     {
00207       double t=exp(-c-exp(logB));
00208       double logt=-c-exp(logB);
00209       double u=t*exp(t);
00210       x0=t*exp(u);
00211       log_x0=logt+u;
00212     }
00213     else if ( (log(.15)<=logB && logB <=log(.35)) ||
00214        ((log(.15)<=logB && logB <=log(.45)) && a>=.3) )
00215     {
00216       double y=-logB;
00217       double v=y-(1-a)*log(y);
00218       x0=y-(1-a)*log(v)-log(1+(1.0-a)/(1.0+v));
00219       log_x0=log(x0);
00220     }
00221     else if (log(.01)<logB && logB < log(.15))
00222     {
00223       double y=-logB;
00224       double v=y-(1-a)*log(y);
00225       x0=y-(1-a)*log(v)-log((v*v+2*(3-a)*v+(2-a)*(3-a))/(v*v +(5-a)*v+2));
00226       log_x0=log(x0);
00227     }
00228     else if (logB < log(.01))
00229     {
00230       double y=-logB;
00231       double v=y-(1-a)*log(y);
00232       x0=y-(1-a)*log(v)-log((v*v+2*(3-a)*v+(2-a)*(3-a))/(v*v +(5-a)*v+2));
00233       log_x0=log(x0);
00234     }
00235     else
00236     {
00237       cerr << "this can't happen" << endl;
00238       ad_exit(1);
00239     }
00240   }
00241   else  if (a>=1.0)
00242   {
00243     const double a0 = 3.31125922108741;
00244     const double b1 = 6.61053765625462;
00245     const double a1 = 11.6616720288968;
00246     const double b2 = 6.40691597760039;
00247     const double a2 = 4.28342155967104;
00248     const double b3 = 1.27364489782223;
00249     const double a3 = .213623493715853;
00250     const double b4 = .03611708101884203;
00251 
00252     int sgn=1;
00253     double logtau;
00254     if (logP< log(0.5))
00255     {
00256       logtau=logP;
00257       sgn=-1;
00258     }
00259     else
00260     {
00261       logtau=logQ;
00262       sgn=1;
00263     }
00264 
00265     double t=sqrt(-2.0*logtau);
00266 
00267     double num = (((a3*t+a2)*t+a1)*t)+a0;
00268     double den = ((((b4*t+b3)*t+b2)*t)+b1)*t+1;
00269     double s=sgn*(t-num/den);
00270     double s2=s*s;
00271     double s3=s2*s;
00272     double s4=s3*s;
00273     double s5=s4*s;
00274     double roota=sqrt(a);
00275     double w=a+s*roota+(s2-1)/3.0+(s3-7.0*s)/(36.*roota)
00276       -(3.0*s4+7.0*s2-16)/(810.0*a)
00277       +(9.0*s5+256.0*s3-433.0*s)/(38880.0*a*roota);
00278     if (logP< log(0.5))
00279     {
00280       if (w>.15*(a+1))
00281       {
00282         x0=w;
00283       }
00284       else
00285       {
00286         double v=logP+gammln(a+1);
00287         double u1=exp(v+w)/a;
00288         double S1=1+u1/(a+1);
00289         double u2=exp((v+u1-log(S1))/a);
00290         double S2=1+u2/(a+1)+u2*u2/((a+1)*(a+2));
00291         double u3=exp((v+u2-log(S2))/a);
00292         double S3=1+u3/(a+1)+u3*u3/((a+1)*(a+2))
00293          + u3*u3*u3/((a+1)*(a+2)*(a+3));
00294         double z=exp((v+u3-log(S3))/a);
00295         if (z<.002*(a+1.0))
00296         {
00297           x0=z;
00298         }
00299         else
00300         {
00301           double sn=Sn(z,a);
00302           double zbar=exp((v+z-log(sn))/a);
00303           x0=zbar*(1.0-(a*log(zbar)-zbar-v+log(sn))/(a-zbar));
00304         }
00305       }
00306       log_x0=log(x0);
00307     }
00308     else
00309     {
00310       double u = -logB +(a-1.0)*log(w)-log(1.0+(1.0-a)/(1+w));
00311       x0=u;
00312       log_x0=log(x0);
00313     }
00314   }
00315   if (a==1.0)
00316   {
00317     x0=-log(1.0-y);
00318     log_x0=log(x0);
00319   }
00320   return log_x0-log(a);
00321 }