ADMB Documentation  11.1.2494
 All Classes Files Functions Variables Typedefs Friends Defines
ccumdbetainv.cpp
Go to the documentation of this file.
00001 
00006 #include <admodel.h>
00007 
00008 static double lnbeta(double a,double b)
00009 {
00010   return gammln(a)+gammln(b)-gammln(a+b);
00011 }
00012 
00013 /*
00014 static int sgn(double z)
00015 {
00016   if (z>=0)
00017     return 1;
00018   else
00019     return -1;
00020 }
00021 */
00022 
00023 double inv_cumd_beta_stable(double a,double b,double y,double eps)
00024 {
00025   double eps1=1.0-eps;
00026 
00027   int icount=0;
00028   if (y<0.0 || y>1.0 || a<=0.0 || b<=0.0 )
00029   {
00030     cerr << "Illegal value in inv_cumd_beta" << endl;
00031     ad_exit(1);
00032   }
00033 
00034   double mu=a/(a+b);
00035   double x=mu;
00036   if (x<=eps)
00037     x=1.1*eps;
00038 
00039   if(x>=eps1)
00040     x=eps1-0.1*eps;
00041 
00042   double s=log(x/(1.0-x));
00043   double lower=log(eps/eps1);
00044   double upper=-lower; // bracket the minimum
00045 
00046   double bet=exp(lnbeta(a,b)); // normalizing constant
00047   double denom=1.0/(betai(a,b,eps1)-betai(a,b,eps));
00048   double finit=betai(a,b,x)-betai(a,b,eps)*denom;
00049 
00050   if (finit>y)
00051   {
00052     upper=s;
00053   }
00054   else if (finit<y)
00055   {
00056     lower=s;
00057   }
00058 
00059   double d=0.0;
00060   do
00061   {
00062     double f,dx; // der of x wrt s
00063     x=1.0/(1.0+exp(-s));  //transform from s to x
00064 
00065     f=(betai(a,b,x)-betai(a,b,eps))*denom;
00066 
00067     dx=exp(-s)/square(1+exp(-s)); // der of x wrt s
00068 
00069     d=y-f;
00070 
00071     if (d<0)   // x is too large
00072     {
00073        if (s<upper)
00074          upper=s;
00075     }
00076     else if (d>0) // x is too small
00077     {
00078       if (s>lower)
00079       {
00080         lower=s;
00081       }
00082     }
00083 
00084     double xa1=pow(x,a-1);
00085     double xb1;
00086     xb1=pow(1-x,b-1);
00087 
00088     double fp=(xa1*xb1/bet)*dx*denom; // derivative of cumulative dist fun
00089                            // wrt x
00090     double h=d/fp;
00091 
00092     double stry=s+h;
00093     if (h<0.0)
00094     {
00095       if (stry<lower)
00096         s=lower+0.5*(s-lower);
00097       else
00098         s=stry;
00099     }
00100     else
00101     if (h>0.0)
00102     {
00103       if (stry>upper)
00104         s=s+0.5*(upper-x);
00105       else
00106         s=stry;
00107     }
00108     icount++;
00109     if (icount>15) break;
00110   }
00111   while(fabs(d)>1.e-12);
00112 
00113   return x;
00114 }
00115 
00116 /*
00117 main()
00118 {
00119   double eps=0.000001;
00120   do
00121   {
00122     double a,b,y;
00123     cin >> a;
00124     cin >> b;
00125     cin >> y;
00126     double x1=inv_cumd_beta(a,b,y);
00127     double x2=inv_cumd_beta_stable(a,b,y,eps);
00128     cout << x1 << " " << x2 << endl;
00129   }
00130   while(1);
00131 }
00132 */