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