ADMB Documentation  11.1.2192
 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()
00163                  +F_xzz * *z.get_u_y()
00164                  +F_x * *z.get_u_yzz());
00165 
00166      double F_zzz=-(F_xxx*cube(*z.get_u_z())+
00167                  +3.0*F_xxz*square(*z.get_u_z())
00168                  +3.0*F_xx* *z.get_u_z() * *z.get_u_zz()
00169                  +3.0*F_xz* *z.get_u_zz()
00170                  +3.0*F_xzz * *z.get_u_z()
00171                  +F_x * *z.get_u_zzz());
00172 
00173      df1b2variable tmp;
00174      double * xd=_y.get_u_dot();
00175      double * yd=_a.get_u_dot();
00176      double * zd=_b.get_u_dot();
00177      double * tmpd=tmp.get_u_dot();
00178      *tmp.get_u()=cx;
00179      for (int i=0;i<df1b2variable::nvar;i++)
00180      {
00181        *tmpd++ = F_x * *xd++ + F_y * *yd++ + F_z * *zd++;
00182      }
00183      if (!df1b2_gradlist::no_derivatives)
00184      {
00185        f1b2gradlist->write_pass1(&_y,&_a,&_b,&tmp,
00186         F_x,F_y,F_z,
00187         F_xx,F_xy,F_xz,F_yy,F_yz,F_zz,F_xxx,F_xxy,F_xxz,F_xyy,F_xyz,
00188         F_xzz,F_yyy,F_yyz,F_yzz,F_zzz);
00189      }
00190 
00191     return tmp;
00192   }
00193 
00205   df3_three_variable betai(const df3_three_variable& a,
00206     const df3_three_variable& b,const  df3_three_variable& x,int maxit)
00207   {
00208     df3_three_variable bt;
00209 
00210     if (value(x) < 0.0 || value(x) > 1.0)
00211       cerr << "Bad x in routine betai" << endl;
00212     if (value(x) == 0.0 || value(x) == 1.0) bt=0.0;
00213     else
00214       bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
00215     if (value(x) < (value(a)+1.0)/(value(a)+value(b)+2.0))
00216       return bt*betacf(a,b,x)/a;
00217     else
00218       return 1.0-bt*betacf(b,a,1.0-x)/b;
00219   }
00220 
00221   df3_three_variable betai(const df3_three_variable& a,
00222     const df3_three_variable& b, double x,int maxit)
00223   {
00224     df3_three_variable bt;
00225 
00226     if (x < 0.0 || x > 1.0) cerr << "Bad x in routine betai" << endl;
00227     if (x == 0.0 || x == 1.0) bt=0.0;
00228     else
00229       bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
00230     if (x < (value(a)+1.0)/(value(a)+value(b)+2.0))
00231       return bt*betacf(a,b,x)/a;
00232     else
00233       return 1.0-bt*betacf(b,a,1.0-x)/b;
00234   }
00235 
00236   df3_three_variable betacf(const df3_three_variable& a,
00237     const  df3_three_variable& b, double x)
00238   {
00239     int m,m2;
00240     df3_three_variable aa,c,d,del,h,qab,qam,qap;
00241 
00242     qab=a+b;
00243     qap=a+1.0;
00244     qam=a-1.0;
00245     c=1.0;
00246     d=1.0-qab*x/qap;
00247     if (fabs(value(d)) < FPMIN) d=FPMIN;
00248     d=1.0/d;
00249     h=d;
00250     for (m=1;m<=MAXIT;m++)
00251     {
00252       m2=2*m;
00253       aa=m*(b-m)*x/((qam+m2)*(a+m2));
00254       d=1.0+aa*d;
00255       if (fabs(value(d)) < FPMIN) d=FPMIN;
00256       c=1.0+aa/c;
00257       if (fabs(value(c)) < FPMIN) c=FPMIN;
00258       d=1.0/d;
00259       h *= d*c;
00260       aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+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 
00267       del=d*c;
00268       h *= del;
00269       if (fabs(value(del)-1.0) < EPS) break;
00270     }
00271     if (m > MAXIT)
00272     {
00273       cerr << "mum interations exceeded " << endl;
00274       ad_exit(1);
00275     }
00276     return h;
00277   }
00278 
00291   df3_three_variable betacf(const df3_three_variable& a,
00292     const df3_three_variable& b, const df3_three_variable& x)
00293   {
00294     int m,m2;
00295     df3_three_variable aa,c,d,del,h,qab,qam,qap;
00296 
00297     qab=a+b;
00298     qap=a+1.0;
00299     qam=a-1.0;
00300     c=1.0;
00301     d=1.0-qab*x/qap;
00302     if (fabs(value(d)) < FPMIN) d=FPMIN;
00303     d=1.0/d;
00304     h=d;
00305     for (m=1;m<=MAXIT;m++)
00306     {
00307       m2=2*m;
00308       aa=m*(b-m)*x/((qam+m2)*(a+m2));
00309       d=1.0+aa*d;
00310       if (fabs(value(d)) < FPMIN) d=FPMIN;
00311       c=1.0+aa/c;
00312       if (fabs(value(c)) < FPMIN) c=FPMIN;
00313       d=1.0/d;
00314       h *= d*c;
00315       aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
00316       d=1.0+aa*d;
00317       if (fabs(value(d)) < FPMIN) d=FPMIN;
00318       c=1.0+aa/c;
00319       if (fabs(value(c)) < FPMIN) c=FPMIN;
00320       d=1.0/d;
00321 
00322       del=d*c;
00323       h *= del;
00324       if (fabs(value(del)-1.0) < EPS) break;
00325     }
00326     if (m > MAXIT)
00327     {
00328       cerr << "mum interations exceeded " << endl;
00329       ad_exit(1);
00330     }
00331     return h;
00332   }
00333 
00334 static df3_three_variable gammlnguts(const df3_three_variable& _z)
00335 {
00336   df3_three_variable x;
00337   const double lpp =0.9189385332046727417803297;
00338   int n=7;
00339   const double c[9]={0.99999999999980993,
00340     676.5203681218851,
00341     -1259.1392167224028,
00342      771.32342877765313,
00343     -176.61502916214059,
00344     12.507343278686905,
00345      -0.13857109526572012,
00346     9.9843695780195716e-6,
00347     1.5056327351493116e-7};
00348   df3_three_variable z=_z-1.0;
00349   x=c[0];
00350   for (int i=1;i<=n+1;i++)
00351   {
00352     df3_three_variable zinv=1.0/(z+i);
00353     x+=c[i]*zinv;
00354   }
00355   df3_three_variable t=z+n+0.5;
00356   df3_three_variable ans= lpp + (z+0.5)*log(t) -t + log(x);
00357   return(ans);
00358 }
00359 
00360 df3_three_variable gammln(const df3_three_variable& z)
00361 {
00362   const double lpi =1.1447298858494001741434272;
00363   const double pi =3.1415926535897932384626432;
00364   if (value(z)<0.5)
00365   {
00366     return lpi - log(sin(pi*z)) - gammlnguts(1.0-z);
00367   }
00368   else
00369   {
00370     return gammlnguts(z);
00371   }
00372 }