ADMB Documentation  11.6rc.3396
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2invcumdbeta.cpp
Go to the documentation of this file.
00001 
00006 #include <df1b2fun.h>
00007 #include <admodel.h>
00008 #include <df3fun.h>
00009 #include "df33fun.h"
00010 #define MAXIT 100
00011 #define EPS 1.0e-9
00012 //#define EPS 3.0e-7
00013 #define FPMIN 1.0e-30
00014   df1b2variable inv_cumd_beta_stable(const df1b2variable& a,
00015     const df1b2variable& b,const df1b2variable& x);
00016   df1b2variable inv_cumd_beta_stable(const df1b2variable& a,
00017     const df1b2variable& b,const df1b2variable& x,double eps);
00018 
00019   df3_three_variable gammln(const df3_three_variable& xx);
00020   df3_three_variable betacf(const df3_three_variable& a,
00021     const df3_three_variable& b,const df3_three_variable& x);
00022   df3_three_variable betacf(const df3_three_variable& a,
00023     const df3_three_variable& b, double x);
00024   df3_three_variable betai(const df3_three_variable& a,
00025     const df3_three_variable& b,const df3_three_variable& x,int maxit);
00026   df3_three_variable betai(const df3_three_variable& a,
00027     const df3_three_variable& b, double x,int maxit);
00028 
00029 /*
00030   static double lnbeta(double a,double b)
00031   {
00032     return gammln(a)+gammln(b)-gammln(a+b);
00033   }
00034 
00035   static int sgn(double z)
00036   {
00037     if (z>=0)
00038       return 1;
00039     else
00040       return -1;
00041   }
00042 */
00043 
00044   double inv_cumd_beta(double _a,double _b,double _y);
00045   double inv_cumd_beta_stable(double _a,double _b,double _y,double eps);
00046 
00047   df1b2variable inv_cumd_beta_stable(const df1b2variable& _a,
00048     const df1b2variable& _b,const df1b2variable& _y,double eps)
00049   {
00050     ADUNCONST(df1b2variable,a);
00051     ADUNCONST(df1b2variable,b);
00052     ADUNCONST(df1b2variable,y);
00053 
00054     double eps1=1.0-eps;
00055     // this gets the inverse without derivatives
00056     double ca=value(a);
00057     double cb=value(b);
00058     double cx=inv_cumd_beta_stable(ca,cb,value(y),eps);
00059 
00060     init_df3_three_variable vx(cx);
00061     init_df3_three_variable va(_a);
00062     init_df3_three_variable vb(_b);
00063 
00064     // this gets the derivatives for the function itself
00065 
00066     df3_three_variable z=(betai(va,vb,vx,25)-betai(va,vb,eps,25))/
00067       (betai(va,vb,eps1,25)-betai(va,vb,eps,25));
00068 
00069     // now solve for the derivatves of the inverse function
00070     double F_x=1.0/(*z.get_u_x());
00071     double F_y=-F_x* *z.get_u_y();
00072     double F_z=-F_x* *z.get_u_z();
00073 
00074     double F_xx=-F_x* *z.get_u_xx()/square(*z.get_u_x());
00075 
00076     double F_xy=-(F_xx * *z.get_u_x() * *z.get_u_y() + F_x * *z.get_u_xy())/
00077                 *z.get_u_x();
00078 
00079     double F_xz=-(F_xx * *z.get_u_x() * *z.get_u_z() + F_x * *z.get_u_xz())/
00080                 *z.get_u_x();
00081     double F_yy=-(F_xx * square(*z.get_u_y())
00082                 +2.0* F_xy* *z.get_u_y()
00083                 +F_x * *z.get_u_yy());
00084 
00085     double F_yz=-( F_xx * *z.get_u_y() * *z.get_u_z()
00086                  + F_x * *z.get_u_yz()
00087                  + F_xy * *z.get_u_z()
00088                  + F_xz * *z.get_u_y());
00089 
00090     double F_zz=-(F_xx * square(*z.get_u_z())
00091                 +2.0* F_xz* *z.get_u_z()
00092                 +F_x * *z.get_u_zz());
00093 
00094     double F_xxx=-(F_x* *z.get_u_xxx()
00095                 +3.0*F_xx* *z.get_u_x() * *z.get_u_xx())
00096                 /cube(*z.get_u_x());
00097     double F_xxy=-(F_xxx * square(*z.get_u_x())* *z.get_u_y()
00098                  + 2.0*F_xx* *z.get_u_x()* *z.get_u_xy()
00099                  + F_xx* *z.get_u_xx() * *z.get_u_y()
00100                  + F_xy * *z.get_u_xx()
00101                  + F_x * *z.get_u_xxy())/
00102                  square(*z.get_u_x());
00103     double F_xxz=-(F_xxx * square(*z.get_u_x())* *z.get_u_z()
00104                  + 2.0*F_xx* *z.get_u_x()* *z.get_u_xz()
00105                  + F_xx* *z.get_u_xx() * *z.get_u_z()
00106                  + F_xz * *z.get_u_xx()
00107                  + F_x * *z.get_u_xxz())/
00108                  square(*z.get_u_x());
00109     double F_xyy=-(F_xxx* *z.get_u_x() *square(*z.get_u_y())
00110                  +2.0* F_xx * *z.get_u_xy() * *z.get_u_y()
00111                  +2.0*F_xxy * *z.get_u_x() * *z.get_u_y()
00112                  + 2.0*F_xy * *z.get_u_xy()
00113                  + F_xx * *z.get_u_x() * *z.get_u_yy()
00114                  + F_x * *z.get_u_xyy())/
00115                  *z.get_u_x();
00116     double F_xyz=-(F_xxx* *z.get_u_x() * *z.get_u_y() * *z.get_u_z()
00117 
00118                  +F_xx * *z.get_u_xy() * *z.get_u_z()
00119                  +F_xx * *z.get_u_xz() * *z.get_u_y()
00120                  +F_xx * *z.get_u_x() * *z.get_u_yz()
00121 
00122                  +F_xxy * *z.get_u_x() * *z.get_u_z()
00123                  +F_xxz * *z.get_u_x() * *z.get_u_y()
00124 
00125                  + F_xy * *z.get_u_xz()
00126                  + F_xz * *z.get_u_xy()
00127 
00128                  + F_x * *z.get_u_xyz())/
00129                  *z.get_u_x();
00130 
00131     double F_xzz=-(F_xxx* *z.get_u_x() *square(*z.get_u_z())
00132                  +2.0* F_xx * *z.get_u_xz() * *z.get_u_z()
00133                  +2.0*F_xxz * *z.get_u_x() * *z.get_u_z()
00134                  + 2.0*F_xz * *z.get_u_xz()
00135                  + F_xx * *z.get_u_x() * *z.get_u_zz()
00136                  + F_x * *z.get_u_xzz())/
00137                  *z.get_u_x();
00138      double F_yyy=-(F_xxx*cube(*z.get_u_y())+
00139                  +3.0*F_xxy*square(*z.get_u_y())
00140                  +3.0*F_xx* *z.get_u_y() * *z.get_u_yy()
00141                  +3.0*F_xy* *z.get_u_yy()
00142                  +3.0*F_xyy * *z.get_u_y()
00143                  +F_x * *z.get_u_yyy());
00144 
00145      double F_yyz=-(F_xxx * square(*z.get_u_y())* *z.get_u_z()
00146                  +2.0*F_xxy * *z.get_u_y() * *z.get_u_z()
00147                  +2.0*F_xyz* *z.get_u_y()
00148                  +F_xxz*square(*z.get_u_y())
00149                  +2.0*F_xx* *z.get_u_y() * *z.get_u_yz()
00150                  +F_xx* *z.get_u_z() * *z.get_u_yy()
00151                  +2.0*F_xy* *z.get_u_yz()
00152                  +F_xz * *z.get_u_yy()
00153                  +F_xyy * *z.get_u_z()
00154                  +F_x * *z.get_u_yyz());
00155      double F_yzz=-(F_xxx * square(*z.get_u_z())* *z.get_u_y()
00156                  +2.0*F_xxz * *z.get_u_z() * *z.get_u_y()
00157                  +2.0*F_xyz* *z.get_u_z()
00158                  +F_xxy*square(*z.get_u_z())
00159                  +2.0*F_xx* *z.get_u_z() * *z.get_u_yz()
00160                  +F_xx* *z.get_u_y() * *z.get_u_zz()
00161                  +2.0*F_xz* *z.get_u_yz()
00162      // +F_xy * *z.get_u_zz()* *z.get_u_y() //** got it!
00163                  +F_xy * *z.get_u_zz() //**
00164                  +F_xzz * *z.get_u_y()
00165                  +F_x * *z.get_u_yzz());
00166 
00167      double F_zzz=-(F_xxx*cube(*z.get_u_z())+
00168                  +3.0*F_xxz*square(*z.get_u_z())
00169                  +3.0*F_xx* *z.get_u_z() * *z.get_u_zz()
00170                  +3.0*F_xz* *z.get_u_zz()
00171                  +3.0*F_xzz * *z.get_u_z()
00172                  +F_x * *z.get_u_zzz());
00173 
00174      df1b2variable tmp;
00175      double * xd=_y.get_u_dot();
00176      double * yd=_a.get_u_dot();
00177      double * zd=_b.get_u_dot();
00178      double * tmpd=tmp.get_u_dot();
00179      *tmp.get_u()=cx;
00180      for (unsigned int i=0;i<df1b2variable::nvar;i++)
00181      {
00182        *tmpd++ = F_x * *xd++ + F_y * *yd++ + F_z * *zd++;
00183      }
00184      if (!df1b2_gradlist::no_derivatives)
00185      {
00186        f1b2gradlist->write_pass1(&_y,&_a,&_b,&tmp,
00187         F_x,F_y,F_z,
00188         F_xx,F_xy,F_xz,F_yy,F_yz,F_zz,
00189   F_xxx,F_xxy,F_xxz,F_xyy,F_xyz,F_xzz,F_yyy,F_yyz,F_yzz,F_zzz);
00190      }
00191      //  cout<<endl<<"-----------"<<endl;
00192      //  cout<<_a<<"\t"<<_b<<"\t"<<_y<<endl;
00193      //  cout<<F_x<<"\t"<<F_y<<"\t"<<F_z<<endl;
00194      //  cout<<F_xx<<"\t"<<F_xy<<"\t"<<F_xz<<"\t"<<F_yy<<"\t"<<F_yz<<"\t"<<F_zz<<endl;
00195      //  cout<<F_xxx<<"\t"<<F_xxy<<"\t"<<F_xxz<<"\t"<<F_xyy<<"\t"<<F_xyz<<"\t"<<F_xzz<<"\t"<<F_yyy<<"\t"<<F_yyz<<"\t"<<F_yzz<<"\t"<<F_zzz<<endl;
00196      //  cout<<endl<<"-----------"<<endl;
00197      //  ad_exit(1);
00198     return tmp;
00199   }
00200 
00212   df3_three_variable betai(const df3_three_variable& a,
00213     const df3_three_variable& b,const  df3_three_variable& x,int maxit)
00214   {
00215     df3_three_variable bt;
00216 
00217     if (value(x) < 0.0 || value(x) > 1.0)
00218       cerr << "Bad x in routine betai" << endl;
00219     if (value(x) == 0.0 || value(x) == 1.0) bt=0.0;
00220     else
00221       bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
00222     if (value(x) < (value(a)+1.0)/(value(a)+value(b)+2.0))
00223       return bt*betacf(a,b,x)/a;
00224     else
00225       return 1.0-bt*betacf(b,a,1.0-x)/b;
00226   }
00227 
00228   df3_three_variable betai(const df3_three_variable& a,
00229     const df3_three_variable& b, double x,int maxit)
00230   {
00231     df3_three_variable bt;
00232 
00233     if (x < 0.0 || x > 1.0) cerr << "Bad x in routine betai" << endl;
00234     if (x == 0.0 || x == 1.0) bt=0.0;
00235     else
00236       bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
00237     if (x < (value(a)+1.0)/(value(a)+value(b)+2.0))
00238       return bt*betacf(a,b,x)/a;
00239     else
00240       return 1.0-bt*betacf(b,a,1.0-x)/b;
00241   }
00242 
00243   df3_three_variable betacf(const df3_three_variable& a,
00244     const  df3_three_variable& b, double x)
00245   {
00246     int m,m2;
00247     df3_three_variable aa,c,d,del,h,qab,qam,qap;
00248 
00249     qab=a+b;
00250     qap=a+1.0;
00251     qam=a-1.0;
00252     c=1.0;
00253     d=1.0-qab*x/qap;
00254     if (fabs(value(d)) < FPMIN) d=FPMIN;
00255     d=1.0/d;
00256     h=d;
00257     for (m=1;m<=MAXIT;m++)
00258     {
00259       m2=2*m;
00260       aa=m*(b-m)*x/((qam+m2)*(a+m2));
00261       d=1.0+aa*d;
00262       if (fabs(value(d)) < FPMIN) d=FPMIN;
00263       c=1.0+aa/c;
00264       if (fabs(value(c)) < FPMIN) c=FPMIN;
00265       d=1.0/d;
00266       h *= d*c;
00267       aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
00268       d=1.0+aa*d;
00269       if (fabs(value(d)) < FPMIN) d=FPMIN;
00270       c=1.0+aa/c;
00271       if (fabs(value(c)) < FPMIN) c=FPMIN;
00272       d=1.0/d;
00273 
00274       del=d*c;
00275       h *= del;
00276       if (fabs(value(del)-1.0) < EPS) break;
00277     }
00278     if (m > MAXIT)
00279     {
00280       cerr << "mum interations exceeded " << endl;
00281       ad_exit(1);
00282     }
00283     return h;
00284   }
00285 
00296   df3_three_variable betacf(const df3_three_variable& a,
00297     const df3_three_variable& b, const df3_three_variable& x)
00298   {
00299     int m,m2;
00300     df3_three_variable aa,c,d,del,h,qab,qam,qap;
00301 
00302     qab=a+b;
00303     qap=a+1.0;
00304     qam=a-1.0;
00305     c=1.0;
00306     d=1.0-qab*x/qap;
00307     if (fabs(value(d)) < FPMIN) d=FPMIN;
00308     d=1.0/d;
00309     h=d;
00310     for (m=1;m<=MAXIT;m++)
00311     {
00312       m2=2*m;
00313       aa=m*(b-m)*x/((qam+m2)*(a+m2));
00314       d=1.0+aa*d;
00315       if (fabs(value(d)) < FPMIN) d=FPMIN;
00316       c=1.0+aa/c;
00317       if (fabs(value(c)) < FPMIN) c=FPMIN;
00318       d=1.0/d;
00319       h *= d*c;
00320       aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
00321       d=1.0+aa*d;
00322       if (fabs(value(d)) < FPMIN) d=FPMIN;
00323       c=1.0+aa/c;
00324       if (fabs(value(c)) < FPMIN) c=FPMIN;
00325       d=1.0/d;
00326 
00327       del=d*c;
00328       h *= del;
00329       if (fabs(value(del)-1.0) < EPS) break;
00330     }
00331     if (m > MAXIT)
00332     {
00333       cerr << "mum interations exceeded " << endl;
00334       ad_exit(1);
00335     }
00336     return h;
00337   }
00338 
00339 static df3_three_variable gammlnguts(const df3_three_variable& _z)
00340 {
00341   df3_three_variable x;
00342   const double lpp =0.9189385332046727417803297;
00343   int n=7;
00344   const double c[9]={0.99999999999980993,
00345     676.5203681218851,
00346     -1259.1392167224028,
00347      771.32342877765313,
00348     -176.61502916214059,
00349     12.507343278686905,
00350      -0.13857109526572012,
00351     9.9843695780195716e-6,
00352     1.5056327351493116e-7};
00353   df3_three_variable z=_z-1.0;
00354   x=c[0];
00355   for (int i=1;i<=n+1;i++)
00356   {
00357     df3_three_variable zinv=1.0/(z+i);
00358     x+=c[i]*zinv;
00359   }
00360   df3_three_variable t=z+n+0.5;
00361   df3_three_variable ans= lpp + (z+0.5)*log(t) -t + log(x);
00362   return(ans);
00363 }
00364 
00365 df3_three_variable gammln(const df3_three_variable& z)
00366 {
00367   const double lpi =1.1447298858494001741434272;
00368   const double pi =3.1415926535897932384626432;
00369   if (value(z)<0.5)
00370   {
00371     return lpi - log(sin(pi*z)) - gammlnguts(1.0-z);
00372   }
00373   else
00374   {
00375     return gammlnguts(z);
00376   }
00377 }
00378 
00379 
00380 df1b2variable qbeta(df1b2variable x, df1b2variable a, df1b2variable b, double eps){
00381   return inv_cumd_beta_stable(a,b,x,eps);
00382 }