ADMB Documentation  11.1.1916
 All Classes Files Functions Variables Typedefs Friends Defines
vgamdev.cpp
Go to the documentation of this file.
00001 /*
00002   $Id: vgamdev.cpp 1708 2014-02-28 21:43:41Z johnoel $
00003 
00004   Author: David Fournier
00005   Copyright (c) 2008, 2009, 2010 Regents of the University of California
00006  */
00007 #include <fvar.hpp>
00008 #define ITMAX 100
00009 //#define EPS 3.0e-7
00010 #define EPS 1.0e-9
00011 #define FPMIN 1.0e-30
00012 static void gcf(double& gammcf,double a,double x,double &gln);
00013 static void gser(double& gamser,double a,double x,double& gln);
00014 
00019   dvariable gamma_deviate(const prevariable& _x,const prevariable& _a)
00020   {
00021     prevariable& x= (prevariable&)(_x);
00022     prevariable& a= (prevariable&)(_a);
00023 
00024     dvariable y=cumd_norm(x);
00025 
00026     y=.9999*y+.00005;
00027 
00028     dvariable z=inv_cumd_gamma(y,a);
00029 
00030     return z;
00031   }
00032 
00033 
00034 static double gammp(double a,double x)
00035 {
00036   double gamser,gammcf,gln;
00037 
00038   if (x < 0.0 || a <= 0.0)
00039     cerr << "Invalid arguments in routine gammp" << endl;
00040   if (x < (a+1.0)) {
00041     gser(gamser,a,x,gln);
00042     return gamser;
00043   } else {
00044     gcf(gammcf,a,x,gln);
00045     return 1.0-gammcf;
00046   }
00047 }
00048 
00055 static void gcf(double& gammcf,double a,double x,double &gln)
00056 {
00057   int i;
00058   double an,b,c,d,del,h;
00059 
00060   gln=gammln(a);
00061   b=x+1.0-a;
00062   c=1.0/FPMIN;
00063   d=1.0/b;
00064   h=d;
00065   for (i=1;i<=ITMAX;i++) {
00066     an = -i*(i-a);
00067     b += 2.0;
00068     d=an*d+b;
00069     if (fabs(d) < FPMIN) d=FPMIN;
00070     c=b+an/c;
00071     if (fabs(c) < FPMIN) c=FPMIN;
00072     d=1.0/d;
00073     del=d*c;
00074     h *= del;
00075     if (fabs(del-1.0) < EPS) break;
00076   }
00077   if (i > ITMAX)
00078     cerr << "a too large, ITMAX too small in gcf" << endl;
00079   gammcf=exp(-x+a*log(x)-(gln))*h;
00080 }
00081 
00088 static void gser(double& gamser,double a,double x,double& gln)
00089 {
00090   int n;
00091   double sum,del,ap;
00092 
00093   gln=gammln(a);
00094   if (x <= 0.0) {
00095     if (x < 0.0)
00096       cerr << "x less than 0 in routine gser" << endl;
00097     gamser=0.0;
00098     return;
00099   } else {
00100     ap=a;
00101     del=sum=1.0/a;
00102     for (n=1;n<=ITMAX;n++) {
00103       ++ap;
00104       del *= x/ap;
00105       sum += del;
00106       if (fabs(del) < fabs(sum)*EPS) {
00107         gamser=sum*exp(-x+a*log(x)-(gln));
00108         return;
00109       }
00110     }
00111     cerr << "a too large, ITMAX too small in routine gser" << endl;
00112     return;
00113   }
00114 }
00115 
00116 static double get_initial_u(double a,double y);
00117 
00118 static double Sn(double x,double a);
00119 
00120 #if defined(USE_LAPLACE)
00121 #include <df1b2fun.h>
00122 df3_two_variable cumd_gamma(const df3_two_variable& x,
00123   const df3_two_variable& a);
00124 
00125 dvariable inv_cumd_gamma(const prevariable& _y,const prevariable& _a)
00126 {
00127   double a=value(_a);
00128   double y=value(_y);
00129   if (a<0.05)
00130   {
00131     cerr << "a musdt be > 0.1" << endl;
00132     ad_exit(1);
00133   }
00134   double u=get_initial_u(a,y);
00135   double h;
00136   int loop_counter=0;
00137   do
00138   {
00139     loop_counter++;
00140     double z=gammp(a,a*exp(u));
00141     double d=y-z;
00142     //cout << d << endl;
00143     double log_fprime=a*log(a)+a*(u-exp(u)) -gammln(a);
00144     double fprime=exp(log_fprime);
00145     h=d/fprime;
00146     u+=h;
00147     if (loop_counter>1000)
00148     {
00149       cerr << "Error in inv_cumd_gamma"
00150         " maximum number of interations exceeded for values"
00151         << endl << "  x = " << y << "  a =  " << a  << "  h =  " << h  << endl;
00152     }
00153   }
00154   while(fabs(h)>1.e-12);
00155 
00156   double x=a*exp(u);
00157 
00158   init_df3_two_variable xx(x);
00159   init_df3_two_variable aa(a);
00160   *xx.get_u_x()=1.0;
00161   *aa.get_u_y()=1.0;
00162 
00163   df3_two_variable z=cumd_gamma(xx,aa);
00164   double F_x=1.0/(*z.get_u_x());
00165   double F_y=-F_x*(*z.get_u_y());
00166 
00167   dvariable vz=0.0;
00168   value(vz)=x;
00169 
00170   gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00171     &(vz.v->x),&(_y.v->x),F_x,&(_a.v->x),F_y);
00172 
00173   return vz;
00174 }
00175 #endif //#if defined(USE_LAPLACE)
00176 
00177 #undef ITMAX
00178 #undef EPS
00179 
00180 static double Sn(double x,double a)
00181 {
00182   int i=1;
00183   double xp=x;
00184   double prod=1.0;
00185   double summ=1.0;
00186   double summand;
00187   do
00188   {
00189     prod*=(a+i);
00190     summand=xp/prod;
00191     if (summand<1.e-4) break;
00192     summ+=summand;
00193     i++;
00194     if (i>50)
00195     {
00196       cerr << "convergence error" << endl;
00197       ad_exit(1);
00198     }
00199   }
00200   while(1);
00201   return summ;
00202 }
00203 
00204 static double get_initial_u(double a,double y)
00205 {
00206   const double c=0.57721;
00207   // note that P = y;
00208   double logP=log(y);
00209   double logQ=log(1-y);
00210   double logB=logQ+gammln(a);
00211   double x0=1.e+100;
00212   double log_x0=1.e+100;
00213 
00214   if (a<1.0)
00215   {
00216     if ( logB>log(.6) || (logB > log(.45) && a>=.3) )
00217     {
00218       double logu;
00219       if (logB+logQ > log(1.e-8))
00220       {
00221         logu=(logP+gammln(1.0+a))/a;
00222       }
00223       else
00224       {
00225         logu=-exp(logQ)/a -c;
00226       }
00227       double u=exp(logu);
00228       x0=u/(1-u/(1.0+a));
00229       double tmp=log(1-u/(1.0+a));
00230       log_x0=logu;
00231       log_x0-=tmp;
00232     }
00233     else if ( a<.3 && log(.35) <= logB && logB <= log(.6) )
00234     {
00235       double t=exp(-c-exp(logB));
00236       double logt=-c-exp(logB);
00237       double u=t*exp(t);
00238       x0=t*exp(u);
00239       log_x0=logt+u;
00240     }
00241     else if ( (log(.15)<=logB && logB <=log(.35)) ||
00242        ((log(.15)<=logB && logB <=log(.45)) && a>=.3) )
00243     {
00244       double y=-logB;
00245       double v=y-(1-a)*log(y);
00246       x0=y-(1-a)*log(v)-log(1+(1.0-a)/(1.0+v));
00247       log_x0=log(x0);
00248     }
00249     else if (log(.01)<logB && logB < log(.15))
00250     {
00251       double y=-logB;
00252       double v=y-(1-a)*log(y);
00253       x0=y-(1-a)*log(v)-log((v*v+2*(3-a)*v+(2-a)*(3-a))/(v*v +(5-a)*v+2));
00254       log_x0=log(x0);
00255     }
00256     else if (logB < log(.01))
00257     {
00258       double y=-logB;
00259       double v=y-(1-a)*log(y);
00260       x0=y-(1-a)*log(v)-log((v*v+2*(3-a)*v+(2-a)*(3-a))/(v*v +(5-a)*v+2));
00261       log_x0=log(x0);
00262     }
00263     else
00264     {
00265       cerr << "this can't happen" << endl;
00266       ad_exit(1);
00267     }
00268   }
00269   else  if (a>=1.0)
00270   {
00271     const double a0 = 3.31125922108741;
00272     const double b1 = 6.61053765625462;
00273     const double a1 = 11.6616720288968;
00274     const double b2 = 6.40691597760039;
00275     const double a2 = 4.28342155967104;
00276     const double b3 = 1.27364489782223;
00277     const double a3 = .213623493715853;
00278     const double b4 = .03611708101884203;
00279 
00280     int sgn=1;
00281     double logtau;
00282     if (logP< log(0.5))
00283     {
00284       logtau=logP;
00285       sgn=-1;
00286     }
00287     else
00288     {
00289       logtau=logQ;
00290       sgn=1;
00291     }
00292 
00293     double t=sqrt(-2.0*logtau);
00294 
00295     double num = (((a3*t+a2)*t+a1)*t)+a0;
00296     double den = ((((b4*t+b3)*t+b2)*t)+b1)*t+1;
00297     double s=sgn*(t-num/den);
00298     double s2=s*s;
00299     double s3=s2*s;
00300     double s4=s3*s;
00301     double s5=s4*s;
00302     double roota=sqrt(a);
00303     double w=a+s*roota+(s2-1)/3.0+(s3-7.0*s)/(36.*roota)
00304       -(3.0*s4+7.0*s2-16)/(810.0*a)
00305       +(9.0*s5+256.0*s3-433.0*s)/(38880.0*a*roota);
00306     if (logP< log(0.5))
00307     {
00308       if (w>.15*(a+1))
00309       {
00310         x0=w;
00311       }
00312       else
00313       {
00314         double v=logP+gammln(a+1);
00315         double u1=exp(v+w)/a;
00316         double S1=1+u1/(a+1);
00317         double u2=exp((v+u1-log(S1))/a);
00318         double S2=1+u2/(a+1)+u2*u2/((a+1)*(a+2));
00319         double u3=exp((v+u2-log(S2))/a);
00320         double S3=1+u3/(a+1)+u3*u3/((a+1)*(a+2))
00321          + u3*u3*u3/((a+1)*(a+2)*(a+3));
00322         double z=exp((v+u3-log(S3))/a);
00323         if (z<.002*(a+1.0))
00324         {
00325           x0=z;
00326         }
00327         else
00328         {
00329           double sn=Sn(z,a);
00330           double zbar=exp((v+z-log(sn))/a);
00331           x0=zbar*(1.0-(a*log(zbar)-zbar-v+log(sn))/(a-zbar));
00332         }
00333       }
00334       log_x0=log(x0);
00335     }
00336     else
00337     {
00338       double u = -logB +(a-1.0)*log(w)-log(1.0+(1.0-a)/(1+w));
00339       x0=u;
00340       log_x0=log(x0);
00341     }
00342   }
00343   if (a==1.0)
00344   {
00345     x0=-log(1.0-y);
00346     log_x0=log(x0);
00347   }
00348   return log_x0-log(a);
00349 }