ADMB Documentation  11.6rc.3396
 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 old_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     cerr<<"y: "<<y<<endl;
00032     cerr<<"a: "<<a<<endl;
00033     cerr<<"b: "<<b<<endl;
00034     ad_exit(1);
00035   }
00036 
00037   double mu=a/(a+b);
00038   double x=mu;
00039   if (x<=eps)
00040     x=1.1*eps;
00041 
00042   if(x>=eps1)
00043     x=eps1-0.1*eps;
00044 
00045   double s=log(x/(1.0-x));
00046   double lower=log(eps/eps1);
00047   double upper=-lower; // bracket the minimum
00048 
00049   double bet=exp(lnbeta(a,b)); // normalizing constant
00050   double denom=1.0/(betai(a,b,eps1)-betai(a,b,eps));
00051   double finit=betai(a,b,x)-betai(a,b,eps)*denom;
00052 
00053   if (finit>y)
00054   {
00055     upper=s;
00056   }
00057   else if (finit<y)
00058   {
00059     lower=s;
00060   }
00061 
00062   double d=0.0;
00063   do
00064   {
00065     double f,dx; // der of x wrt s
00066     x=1.0/(1.0+exp(-s));  //transform from s to x
00067 
00068     f=(betai(a,b,x)-betai(a,b,eps))*denom;
00069 
00070     dx=exp(-s)/square(1+exp(-s)); // der of x wrt s
00071 
00072     d=y-f;
00073 
00074     if (d<0)   // x is too large
00075     {
00076        if (s<upper)
00077          upper=s;
00078     }
00079     else if (d>0) // x is too small
00080     {
00081       if (s>lower)
00082       {
00083         lower=s;
00084       }
00085     }
00086 
00087     double xa1=pow(x,a-1);
00088     double xb1;
00089     xb1=pow(1-x,b-1);
00090 
00091     double fp=(xa1*xb1/bet)*dx*denom; // derivative of cumulative dist fun
00092                            // wrt x
00093     double h=d/fp;
00094 
00095     double stry=s+h;
00096     if (h<0.0)
00097     {
00098       if (stry<lower)
00099         s=lower+0.5*(s-lower);
00100       else
00101         s=stry;
00102     }
00103     else
00104     if (h>0.0)
00105     {
00106       if (stry>upper)
00107         s=s+0.5*(upper-x);
00108       else
00109         s=stry;
00110     }
00111     icount++;
00112     if (icount>15) break;
00113   }
00114   while(fabs(d)>1.e-12);
00115 
00116   return x;
00117 }
00118 
00119 
00120 double inv_cumd_beta_stable(double a,double b,double y,double eps)
00121 {
00122   //double eps1=1.0-eps;
00123 
00124   int icount=0;
00125   if (y<0.0 || y>1.0 || a<=0.0 || b<=0.0 )
00126   {
00127     cerr << "Illegal value in inv_cumd_beta" << endl;
00128     cerr<<"y: "<<y<<endl;
00129     cerr<<"a: "<<a<<endl;
00130     cerr<<"b: "<<b<<endl;
00131     ad_exit(1);
00132   }
00133 
00134   double mu=a/(a+b);
00135   double x=mu;
00136 
00137   double lower=0.0;
00138   double upper=1.0; // bracket the minimum
00139 
00140   double bet=exp(lnbeta(a,b)); // normalizing constant
00141   double finit=betai(a,b,x);
00142 
00143   if (finit>y)
00144   {
00145     upper=x;
00146   }
00147   else if (finit<y)
00148   {
00149     lower=x;
00150   }
00151 
00152   double d=0.0;
00153   do
00154   {
00155     double f; // der of x wrt s
00156     f=betai(a,b,x);
00157     d=y-f;
00158     if (d<0)   // x is too large
00159     {
00160        if (x<upper)
00161          upper=x;
00162     }
00163     else if (d>0) // x is too small
00164     {
00165       if (x>lower)
00166       {
00167         lower=x;
00168       }
00169     }
00170     double xa1=pow(x,a-1);
00171     double xb1=pow(1-x,b-1);
00172 
00173     double fp=(xa1*xb1/bet); // derivative of cumulative dist fun wrt x
00174      
00175     double h=d/fp;
00176     double stry=x+h;
00177     if (h<0.0)
00178     {
00179       if (stry<lower)
00180         x=lower+0.5*(x-lower);
00181       else
00182         x=stry;
00183     }
00184     else
00185     if (h>0.0)
00186     {
00187       if (stry>upper)
00188         x=x+0.5*(upper-x);
00189       else
00190         x=stry;
00191     }
00192     icount++;
00193     if (icount>50) break;
00194   }
00195   while(fabs(d)>1.e-16);
00196   return x;
00197 }
00198 
00199 
00200 double qbeta(double x, double a, double b, double eps){
00201   return inv_cumd_beta_stable(a,b,x,eps);
00202 }
00203 
00204 
00205 
00206 /*
00207 main()
00208 {
00209   double eps=0.000001;
00210   do
00211   {
00212     double a,b,y;
00213     cin >> a;
00214     cin >> b;
00215     cin >> y;
00216     double x1=inv_cumd_beta(a,b,y);
00217     double x2=inv_cumd_beta_stable(a,b,y,eps);
00218     cout << x1 << " " << x2 << endl;
00219   }
00220   while(1);
00221 }
00222 */