ccumdbetainv_incbet.cpp
Go to the documentation of this file.
```00001 /*
00002  * \$Id: ccumdbetainv_incbet.cpp 422 2012-04-12 23:02:36Z johnoel \$
00003  * Author: Unknown
00004  */
00010
00015 static double lnbeta(double a,double b)
00016 {
00017   return gammln(a)+gammln(b)-gammln(a+b);
00018 }
00019
00024 double inv_cumd_beta_stable(double a,double b,double y,double eps)
00025 {
00026   double eps1=1.0-eps;
00027
00028   int icount=0;
00029   if (y<0.0 || y>1.0 || a<=0.0 || b<=0.0 )
00030   {
00031     cerr << "Illegal value in inv_cumd_beta" << endl;
00033   }
00034
00035   double mu=a/(a+b);
00036   double x=mu;
00037   if (x<=eps)
00038     x=1.1*eps;
00039
00040   if(x>=eps1)
00041     x=eps1-0.1*eps;
00042
00043   double s=log(x/(1.0-x));
00044   double lower=log(eps/eps1);
00045   double upper=-lower; // bracket the minimum
00046
00047   double bet=exp(lnbeta(a,b)); // normalizing constant
00048   double denom=1.0/(incbet(a,b,eps1)-incbet(a,b,eps));
00049   double finit=incbet(a,b,x)-incbet(a,b,eps)*denom;
00050
00051   if (finit>y)
00052   {
00053     upper=s;
00054   }
00055   else if (finit<y)
00056   {
00057     lower=s;
00058   }
00059
00060   double d=0.0;
00061   do
00062   {
00063     double f,dx; // der of x wrt s
00064     x=1.0/(1.0+exp(-s));  //transform from s to x
00065
00066     f=(incbet(a,b,x)-incbet(a,b,eps))*denom;
00067
00068     dx=exp(-s)/square(1+exp(-s)); // der of x wrt s
00069
00070     d=y-f;
00071
00072     if (d<0)   // x is too large
00073     {
00074        if (s<upper)
00075          upper=s;
00076     }
00077     else if (d>0) // x is too small
00078     {
00079       if (s>lower)
00080       {
00081         lower=s;
00082       }
00083     }
00084
00085     double xa1=pow(x,a-1);
00086     double xb1;
00087     xb1=pow(1-x,b-1);
00088
00089     double fp=(xa1*xb1/bet)*dx*denom; // derivative of cumulative dist fun
00090                            // wrt x
00091     double h=d/fp;
00092
00093     double stry=s+h;
00094     if (h<0.0)
00095     {
00096       if (stry<lower)
00097         s=lower+0.5*(s-lower);
00098       else
00099         s=stry;
00100     }
00101     else
00102     if (h>0.0)
00103     {
00104       if (stry>upper)
00105         s=s+0.5*(upper-x);
00106       else
00107         s=stry;
00108     }
00109     icount++;
00110     if (icount>15) break;
00111   }
00112   while(fabs(d)>1.e-12);
00113
00114   return x;
00115 }
00116
```