ADMB Documentation  11.1x.2728
 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 
00073 void gcf(const df3_two_variable& _gammcf,const df3_two_variable& a,
00074   const df3_two_variable& x,const df3_two_variable& _gln)
00075 {
00076   ADUNCONST(df3_two_variable,gln)
00077   ADUNCONST(df3_two_variable,gammcf)
00078   int i;
00079   df3_two_variable an,b,c,d,del,h;
00080 
00081   gln=gammln(a);
00082   b=x+1.0-a;
00083   c=1.0/FPMIN;
00084   d=1.0/b;
00085   h=d;
00086   for (i=1;i<=ITMAX;i++) {
00087     an = -i*(i-a);
00088     b += 2.0;
00089     d=an*d+b;
00090     if (fabs(value(d)) < FPMIN) d=FPMIN;
00091     c=b+an/c;
00092     if (fabs(value(c)) < FPMIN) c=FPMIN;
00093     d=1.0/d;
00094     del=d*c;
00095     h *= del;
00096     if (fabs(value(del)-1.0) < EPS) break;
00097   }
00098   if (i > ITMAX)
00099     cerr << "a too large, ITMAX too small in gcf" << endl;
00100   gammcf=exp(-x+a*log(x)-(gln))*h;
00101 }
00102 
00111 void gser(const df3_two_variable& _gamser,const df3_two_variable& a,
00112   const df3_two_variable& x,const df3_two_variable& _gln)
00113 {
00114   int n;
00115   ADUNCONST(df3_two_variable,gln)
00116   ADUNCONST(df3_two_variable,gamser)
00117   df3_two_variable sum,del,ap;
00118 
00119   gln=gammln(a);
00120 
00121   if (value(x) <= 0.0) {
00122     if (value(x) < 0.0)
00123       cerr << "x less than 0 in routine gser" << endl;
00124     gamser=0.0;
00125     return;
00126   }
00127   else
00128   {
00129     ap=a;
00130     del=sum=1.0/a;
00131     for (n=1;n<=ITMAX;n++) {
00132       ap+=1.0;
00133       del *= x/ap;
00134       sum += del;
00135       if (fabs(value(del)) < fabs(value(sum))*EPS) {
00136         gamser=sum*exp(-x+a*log(x)-(gln));
00137         return;
00138       }
00139     }
00140     cerr << "a too large, ITMAX too small in routine gser" << endl;
00141     return;
00142   }
00143 }
00144 
00145 
00146 df3_two_variable cumd_gamma(const df3_two_variable& x,
00147   const df3_two_variable& a)
00148 {
00149   df3_two_variable gamser,gammcf,gln;
00150 
00151   if (value(x) < 0.0 || value(a) <= 0.0)
00152     cerr << "Invalid arguments in routine gammp" << endl;
00153   if (value(x) < (value(a)+1.0)) {
00154     gser(gamser,a,x,gln);
00155     return gamser;
00156   } else {
00157     gcf(gammcf,a,x,gln);
00158     return 1.0-gammcf;
00159   }
00160 }
00161 df3_two_variable cumd_exponential(const df3_two_variable& x,
00162   const df3_two_variable& a)
00163 {
00164   df3_two_variable tmp;
00165   if (value(x)<=0)
00166     return 0.5*exp(x);
00167   else
00168     return 1.0-0.5*exp(-x);
00169 }
00170 df3_two_variable cumd_cauchy(const df3_two_variable& x,
00171   const df3_two_variable& a)
00172 {
00173   return atan(x/a);
00174 }