ADMB Documentation  11.6rc.3405
 All Classes Files Functions Variables Typedefs Friends Defines
df3gammp.cpp
Go to the documentation of this file.
00001 
00007 #include <df1b2fun.h>
00008 #define ITMAX 100
00009 #define EPS 1.0e-9
00010 //#define EPS 3.0e-7
00011 #define FPMIN 1.0e-30
00012 double get_values(double x,double y,int print_switch);
00013 
00014 
00015 df1b2variable log_negbinomial_density(double x,const df1b2variable& _xmu,
00016   const df1b2variable& _xtau)
00017 {
00018   ADUNCONST(df1b2variable,xmu)
00019   ADUNCONST(df1b2variable,xtau)
00020   init_df3_two_variable mu(xmu);
00021   init_df3_two_variable tau(xtau);
00022   *mu.get_u_x()=1.0;
00023   *tau.get_u_y()=1.0;
00024   if (value(tau)-1.0<0.0)
00025   {
00026     cerr << "tau <=1 in log_negbinomial_density " << endl;
00027     ad_exit(1);
00028   }
00029   df3_two_variable r=mu/(1.e-120+(tau-1.0));
00030   df3_two_variable tmp;
00031   tmp=gammln(x+r)-gammln(r) -gammln(x+1)
00032     +r*log(r)+x*log(mu)-(r+x)*log(r+mu);
00033   df1b2variable tmp1;
00034   tmp1=tmp;
00035   return tmp1;
00036 }
00037 
00046 df3_two_variable gammln(const df3_two_variable& xx)
00047 {
00048   df3_two_variable x,tmp,ser,tmp1;
00049   static double cof[6]={76.18009173,-86.50532033,24.01409822,
00050     -1.231739516,0.120858003e-2,-0.536382e-5};
00051   int j;
00052   x=xx-1.0;
00053   tmp=x+5.5;
00054   tmp -= (x+0.5)*log(tmp);
00055   ser=1.0;
00056   for (j=0;j<=5;j++)
00057   {
00058     x += 1.0;
00059     ser += cof[j]/x;
00060   }
00061   return -tmp+log(2.50662827465*ser);
00062 }
00063 
00064 
00071 void gcf(const df3_two_variable& _gammcf,const df3_two_variable& a,
00072   const df3_two_variable& x,const df3_two_variable& _gln)
00073 {
00074   ADUNCONST(df3_two_variable,gln)
00075   ADUNCONST(df3_two_variable,gammcf)
00076   int i;
00077   df3_two_variable an,b,c,d,del,h;
00078 
00079   gln=gammln(a);
00080   b=x+1.0-a;
00081   c=1.0/FPMIN;
00082   d=1.0/b;
00083   h=d;
00084   for (i=1;i<=ITMAX;i++) {
00085     an = -i*(i-a);
00086     b += 2.0;
00087     d=an*d+b;
00088     if (fabs(value(d)) < FPMIN) d=FPMIN;
00089     c=b+an/c;
00090     if (fabs(value(c)) < FPMIN) c=FPMIN;
00091     d=1.0/d;
00092     del=d*c;
00093     h *= del;
00094     if (fabs(value(del)-1.0) < EPS) break;
00095   }
00096   if (i > ITMAX)
00097     cerr << "a too large, ITMAX too small in gcf" << endl;
00098   gammcf=exp(-x+a*log(x)-(gln))*h;
00099 }
00100 
00107 void gser(const df3_two_variable& _gamser,const df3_two_variable& a,
00108   const df3_two_variable& x,const df3_two_variable& _gln)
00109 {
00110   int n;
00111   ADUNCONST(df3_two_variable,gln)
00112   ADUNCONST(df3_two_variable,gamser)
00113   df3_two_variable sum,del,ap;
00114 
00115   gln=gammln(a);
00116 
00117   if (value(x) <= 0.0) {
00118     if (value(x) < 0.0)
00119       cerr << "x less than 0 in routine gser" << endl;
00120     gamser=0.0;
00121     return;
00122   }
00123   else
00124   {
00125     ap=a;
00126     del=sum=1.0/a;
00127     for (n=1;n<=ITMAX;n++) {
00128       ap+=1.0;
00129       del *= x/ap;
00130       sum += del;
00131       if (fabs(value(del)) < fabs(value(sum))*EPS) {
00132         gamser=sum*exp(-x+a*log(x)-(gln));
00133         return;
00134       }
00135     }
00136     cerr << "a too large, ITMAX too small in routine gser" << endl;
00137     return;
00138   }
00139 }
00140 
00141 
00142 df3_two_variable cumd_gamma(const df3_two_variable& x,
00143   const df3_two_variable& a)
00144 {
00145   df3_two_variable gamser,gammcf,gln;
00146 
00147   if (value(x) < 0.0 || value(a) <= 0.0)
00148     cerr << "Invalid arguments in routine gammp" << endl;
00149   if (value(x) < (value(a)+1.0)) {
00150     gser(gamser,a,x,gln);
00151     return gamser;
00152   } else {
00153     gcf(gammcf,a,x,gln);
00154     return 1.0-gammcf;
00155   }
00156 }
00157 df3_two_variable cumd_exponential(const df3_two_variable& x,
00158   const df3_two_variable& a)
00159 {
00160   df3_two_variable tmp;
00161   if (value(x)<=0)
00162     return 0.5*exp(x);
00163   else
00164     return 1.0-0.5*exp(-x);
00165 }
00166 df3_two_variable cumd_cauchy(const df3_two_variable& x,
00167   const df3_two_variable& a)
00168 {
00169   return atan(x/a);
00170 }