ADMB Documentation  11.1.1932
 All Classes Files Functions Variables Typedefs Friends Defines
vcumdbetainv.cpp
Go to the documentation of this file.
00001 
00007 #if defined(USE_LAPLACE)
00008 #  include <df1b2fun.h>
00009 #endif
00010 
00011 #include <admodel.h>
00012 #include "df12fun.h"
00013 
00014 //#define ADUNCONST(type,obj) type & obj = (type&) _##obj;
00015 
00016 static double lnbeta(double a,double b)
00017 {
00018   return gammln(a)+gammln(b)-gammln(a+b);
00019 }
00020 
00021 /*
00022 static int sgn(double z)
00023 {
00024   if (z>=0)
00025     return 1;
00026   else
00027     return -1;
00028 }
00029 */
00030 
00031 df1_two_variable betai(const df1_two_variable& a,
00032   const df1_two_variable& b,double x,int maxit=100);
00033 
00034 
00035 dvariable inv_cumd_beta_stable(const prevariable& _a,const prevariable& _b,
00036   const prevariable& _y,double eps)
00037 {
00038   ADUNCONST(prevariable,a);
00039   ADUNCONST(prevariable,b);
00040   ADUNCONST(prevariable,y);
00041 
00042   double eps1=1.0-eps;
00043   // this gets the inverse without derivatives
00044   double ca=value(a);
00045   double cb=value(b);
00046   double cx=inv_cumd_beta_stable(ca,cb,value(y),eps);
00047 
00048   init_df1_two_variable va(_a);
00049   init_df1_two_variable vb(_b);
00050 
00051   // this gets the derivatives for the function itself
00052 
00053   df1_two_variable z=(betai(va,vb,cx)-betai(va,vb,eps))/
00054     (betai(va,vb,eps1)-betai(va,vb,eps));
00055 
00056   double dga=*z.get_u_x();
00057   double dgb=*z.get_u_y();
00058 
00059   double denom=1.0/(betai(ca,cb,eps1)-betai(ca,cb,eps));
00060   double dgx=pow(cx,value(a)-1.0)*pow(1.0-cx,value(b)-1.0)/
00061     exp(lnbeta(value(a),value(b)))*denom;
00062 
00063   // now solve for the derivatves of the inverse function
00064 
00065   double dfx=1.0/dgx;
00066   double dfa=-dfx*dga;
00067   double dfb=-dfx*dgb;
00068 
00069   dvariable tmp;
00070   value(tmp)=cx;
00071 
00072   gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation3ind,
00073     &(value(tmp)) ,&(value(_a)),dfa ,&(value(_b)),dfb ,&(value(_y)),dfx);
00074 
00075   return tmp;
00076 }
00077 
00078 df1_two_variable betacf(const df1_two_variable& a,const df1_two_variable& b,
00079   double x);
00080 
00081 
00082 df1_two_variable betai(const df1_two_variable& a,
00083   const df1_two_variable& b,double x,int maxit)
00084 {
00085   df1_two_variable bt;
00086 
00087   if (x < 0.0 || x > 1.0) cerr << "Bad x in routine betai" << endl;
00088   double z=1.0-x;
00089   if (x == 0.0 || x==1.0)
00090     bt=0.0;
00091   else
00092   {
00093     bt=exp(gammln(a+b) -gammln(a) -gammln(b) +a*log(x) +b*log(z));
00094   }
00095   if (x < (value(a)+1.0)/(value(a)+value(b)+2.0))
00096     return bt*betacf(a,b,x)/a;
00097   else
00098     return 1.0-bt*betacf(b,a,1.0-x)/b;
00099 }
00100 
00101 /*
00102 main()
00103 {
00104   double a,b,x;
00105 
00106   a=2.0;
00107   b=3.0;
00108   x=.80;
00109   do
00110   {
00111     cin >> a;
00112     cin >> b;
00113     cin >> x;
00114     if (x<0 ) exit(1);
00115     independent_variables xx(1,3);
00116     dvector g(1,3);
00117     dvector gg(1,3);
00118     xx(1)=a;
00119     xx(2)=b;
00120     xx(3)=x;
00121     gradient_structure gs(10000);
00122     {
00123       dvar_vector vx=dvar_vector(xx);
00124       dvariable y=inv_cumd_beta(vx(1),vx(2),vx(3));
00125       cout << y << endl;
00126       gradcalc(3,g);
00127       double h=1.e-6;
00128       gg(1)=(inv_cumd_beta(a+h,b,x)-inv_cumd_beta(a-h,b,x))/(2.0*h);
00129       gg(2)=(inv_cumd_beta(a,b+h,x)-inv_cumd_beta(a,b-h,x))/(2.0*h);
00130       gg(3)=(inv_cumd_beta(a,b,x+h)-inv_cumd_beta(a,b,x-h))/(2.0*h);
00131       cout << g << endl  << gg << " " << endl << g-gg << endl;
00132     }
00133     {
00134       dvar_vector vx=dvar_vector(xx);
00135       dvariable y=inv_cumd_beta_stable(vx(1),vx(2),vx(3),1.e-7);
00136       cout << y << endl;
00137       gradcalc(3,g);
00138       double h=1.e-6;
00139       gg(1)=(inv_cumd_beta_stable(a+h,b,x,1.e-7)-inv_cumd_beta_stable(a-h,b,x,1.e-7))/
00140         (2.0*h);
00141       gg(2)=(inv_cumd_beta_stable(a,b+h,x,1.e-7)-inv_cumd_beta_stable(a,b-h,x,1.e-7))/
00142         (2.0*h);
00143       gg(3)=(inv_cumd_beta_stable(a,b,x+h,1.e-7)-inv_cumd_beta_stable(a,b,x-h,1.e-7))/
00144         (2.0*h);
00145       cout << g << endl  << gg << " " << endl << g-gg << endl;
00146     }
00147   }
00148   while(1);
00149 }
00150 */
00151 
00152 const double tiny=1.0e-8;
00153 const double maxn=150;
00154 const double lowerbd=1.0e-40;
00155 
00156 df1_two_variable betacf(const df1_two_variable& a, const df1_two_variable& b,
00157   double x)
00158 {
00159   df1_two_variable v;
00160   df1_two_variable aa;
00161   df1_two_variable h;
00162   df1_two_variable up;
00163   df1_two_variable yy;
00164   df1_two_variable um;
00165   df1_two_variable d;
00166   df1_two_variable ssum;
00167 
00168   up=a+1.0;
00169   um=a-1.0;
00170   v=1.0;
00171   ssum=a+b;
00172   d=1.0-ssum*x/up;
00173   if (fabs(value(d)) < lowerbd) d=lowerbd;
00174   d=1.0/d;
00175   h=d;
00176   int m;
00177   for (m=1;m<=maxn;m++)
00178   {
00179     int m2=2*m;
00180     aa=m*(b-m)*x/((um+m2)*(a+m2));
00181     d=1.0+aa*d;
00182     if (fabs(value(d)) < lowerbd) d=lowerbd;
00183     v=1.0+aa/v;
00184     if (fabs(value(v)) < lowerbd) v=lowerbd;
00185     d=1.0/d;
00186     h *= d*v;
00187     aa = -(a+m)*(ssum+m)*x/((a+m2)*(up+m2));
00188     d=1.0+aa*d;
00189     if (fabs(value(d)) < lowerbd) d=lowerbd;
00190     v=1.0+aa/v;
00191     if (fabs(value(v)) < lowerbd) v=lowerbd;
00192     d=1.0/d;
00193 
00194     yy=d*v;
00195     h *= yy;
00196     if (fabs(value(yy)-1.0) < tiny) break;
00197   }
00198   if (m > maxn)
00199   {
00200     cerr << "num interations exceeded " << endl;
00201     ad_exit(1);
00202   }
00203   return h;
00204 }
00205 
00206 /*
00207 
00208 df1_two_variable gammln(const df1_two_variable& xx)
00209 {
00210   df1_two_variable x,tmp,ser;
00211   static double cof[6]={76.18009173,-86.50532033,24.01409822,
00212     -1.231739516,0.120858003e-2,-0.536382e-5};
00213   int j;
00214   x=xx-1.0;
00215   tmp=x+5.5;
00216   tmp -= (x+0.5)*log(tmp);
00217   ser=1.0;
00218   for (j=0;j<=5;j++)
00219   {
00220     x += 1.0;
00221     ser += cof[j]/x;
00222   }
00223   return -tmp+log(2.50662827465*ser);
00224 }
00225 */
00226 
00227 
00228 static df1_two_variable gammlnguts(const df1_two_variable& _z)
00229 {
00230   df1_two_variable x;
00231   const double lpp =0.9189385332046727417803297;
00232   int n=7;
00233   const double c[9]={0.99999999999980993,
00234     676.5203681218851,
00235     -1259.1392167224028,
00236      771.32342877765313,
00237     -176.61502916214059,
00238     12.507343278686905,
00239      -0.13857109526572012,
00240     9.9843695780195716e-6,
00241     1.5056327351493116e-7};
00242   df1_two_variable z=_z-1.0;
00243   x=c[0];
00244   for (int i=1;i<=n+1;i++)
00245   {
00246     df1_two_variable zinv=1.0/(z+i);
00247     x+=c[i]*zinv;
00248   }
00249   df1_two_variable t=z+n+0.5;
00250   df1_two_variable ans= lpp + (z+0.5)*log(t) -t + log(x);
00251   return(ans);
00252 }
00253 
00254 df1_two_variable gammln(const df1_two_variable& z)
00255 {
00256   const double lpi =1.1447298858494001741434272;
00257   const double pi =3.1415926535897932384626432;
00258   if (value(z)<0.5)
00259   {
00260     return lpi - log(sin(pi*z)) - gammlnguts(1.0-z);
00261   }
00262   else
00263   {
00264     return gammlnguts(z);
00265   }
00266 }
00267