ADMB Documentation  11.1.2547
 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 get_initial_u(double a,double y)
00153 {
00154   const double c=0.57721;
00155   // note that P = y;
00156   double logP=log(y);
00157   double logQ=log(1-y);
00158   double logB=logQ+gammln(a);
00159   double x0=1.e+100;
00160   double log_x0=1.e+100;
00161 
00162   if (a<1.0)
00163   {
00164     if ( logB>log(.6) || (logB > log(.45) && a>=.3) )
00165     {
00166       double logu;
00167       if (logB+logQ > log(1.e-8))
00168       {
00169         logu=(logP+gammln(1.0+a))/a;
00170       }
00171       else
00172       {
00173         logu=-exp(logQ)/a -c;
00174       }
00175       double u=exp(logu);
00176       x0=u/(1-u/(1.0+a));
00177       double tmp=log(1-u/(1.0+a));
00178       log_x0=logu;
00179       log_x0-=tmp;
00180     }
00181     else if ( a<.3 && log(.35) <= logB && logB <= log(.6) )
00182     {
00183       double t=exp(-c-exp(logB));
00184       double logt=-c-exp(logB);
00185       double u=t*exp(t);
00186       x0=t*exp(u);
00187       log_x0=logt+u;
00188     }
00189     else if ( (log(.15)<=logB && logB <=log(.35)) ||
00190        ((log(.15)<=logB && logB <=log(.45)) && a>=.3) )
00191     {
00192       double y=-logB;
00193       double v=y-(1-a)*log(y);
00194       x0=y-(1-a)*log(v)-log(1+(1.0-a)/(1.0+v));
00195       log_x0=log(x0);
00196     }
00197     else if (log(.01)<logB && logB < log(.15))
00198     {
00199       double y=-logB;
00200       double v=y-(1-a)*log(y);
00201       x0=y-(1-a)*log(v)-log((v*v+2*(3-a)*v+(2-a)*(3-a))/(v*v +(5-a)*v+2));
00202       log_x0=log(x0);
00203     }
00204     else if (logB < log(.01))
00205     {
00206       double y=-logB;
00207       double v=y-(1-a)*log(y);
00208       x0=y-(1-a)*log(v)-log((v*v+2*(3-a)*v+(2-a)*(3-a))/(v*v +(5-a)*v+2));
00209       log_x0=log(x0);
00210     }
00211     else
00212     {
00213       cerr << "this can't happen" << endl;
00214       ad_exit(1);
00215     }
00216   }
00217   else  if (a>=1.0)
00218   {
00219     const double a0 = 3.31125922108741;
00220     const double b1 = 6.61053765625462;
00221     const double a1 = 11.6616720288968;
00222     const double b2 = 6.40691597760039;
00223     const double a2 = 4.28342155967104;
00224     const double b3 = 1.27364489782223;
00225     const double a3 = .213623493715853;
00226     const double b4 = .03611708101884203;
00227 
00228     int sgn=1;
00229     double logtau;
00230     if (logP< log(0.5))
00231     {
00232       logtau=logP;
00233       sgn=-1;
00234     }
00235     else
00236     {
00237       logtau=logQ;
00238       sgn=1;
00239     }
00240 
00241     double t=sqrt(-2.0*logtau);
00242 
00243     double num = (((a3*t+a2)*t+a1)*t)+a0;
00244     double den = ((((b4*t+b3)*t+b2)*t)+b1)*t+1;
00245     double s=sgn*(t-num/den);
00246     double s2=s*s;
00247     double s3=s2*s;
00248     double s4=s3*s;
00249     double s5=s4*s;
00250     double roota=sqrt(a);
00251     double w=a+s*roota+(s2-1)/3.0+(s3-7.0*s)/(36.*roota)
00252       -(3.0*s4+7.0*s2-16)/(810.0*a)
00253       +(9.0*s5+256.0*s3-433.0*s)/(38880.0*a*roota);
00254     if (logP< log(0.5))
00255     {
00256       if (w>.15*(a+1))
00257       {
00258         x0=w;
00259       }
00260       else
00261       {
00262         double v=logP+gammln(a+1);
00263         double u1=exp(v+w)/a;
00264         double S1=1+u1/(a+1);
00265         double u2=exp((v+u1-log(S1))/a);
00266         double S2=1+u2/(a+1)+u2*u2/((a+1)*(a+2));
00267         double u3=exp((v+u2-log(S2))/a);
00268         double S3=1+u3/(a+1)+u3*u3/((a+1)*(a+2))
00269          + u3*u3*u3/((a+1)*(a+2)*(a+3));
00270         double z=exp((v+u3-log(S3))/a);
00271         if (z<.002*(a+1.0))
00272         {
00273           x0=z;
00274         }
00275         else
00276         {
00277           double sn=Sn(z,a);
00278           double zbar=exp((v+z-log(sn))/a);
00279           x0=zbar*(1.0-(a*log(zbar)-zbar-v+log(sn))/(a-zbar));
00280         }
00281       }
00282       log_x0=log(x0);
00283     }
00284     else
00285     {
00286       double u = -logB +(a-1.0)*log(w)-log(1.0+(1.0-a)/(1+w));
00287       x0=u;
00288       log_x0=log(x0);
00289     }
00290   }
00291   if (a==1.0)
00292   {
00293     x0=-log(1.0-y);
00294     log_x0=log(x0);
00295   }
00296   return log_x0-log(a);
00297 }