ADMB Documentation  11.1.2491
 All Classes Files Functions Variables Typedefs Friends Defines
cumdist.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: cumdist.cpp 1112 2013-07-12 21:41:41Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <fvar.hpp>
00012 
00013 double inv_cumd_norm(const double& x);
00014 double cumd_norm(const double& x);
00015 
00020 double normal_tail_right(const double& x)
00021 {
00022   const double a3=5;
00023   const double a4=9;
00024   const double a5=129;
00025   double x2=x*x;
00026   double z=0.3989422804*exp(-.5*x2);
00027   double b1=x2+2;
00028   double b2=b1*(x2+4);
00029   double b3=b2*(x2+6);
00030   double b4=b3*(x2+8);
00031   double b5=b4*(x2+10);
00032   double tmp=z/x*(1.0 -1.0/b1 +1.0/b2 - a3/b3 +a4/b4 -a5/b5);
00033   return tmp;
00034 }
00035 
00040 double inv_cumd_norm_inner(double x)
00041 {
00042   const double c0=2.515517;
00043   const double c1=0.802853;
00044   const double c2=0.010328;
00045   const double d1=1.432788;
00046   const double d2=0.189269;
00047   const double d3=0.001308;
00048   if (x<=0 || x>=1.0)
00049   {
00050     cerr << "Illegal argument to inv_cumd_norm = " << x << endl;
00051     return 0;
00052   }
00053 
00054   if (x<0.5)
00055   {
00056     double t = sqrt(-2.*log(x));
00057     double p=((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1)-t;
00058     return p;
00059   }
00060   else if (x==0.5)
00061   {
00062     return 0.0;
00063   }
00064   else
00065   {
00066     double y=1.-x;
00067     //double t = sqrt(-log(y*y));
00068     double t = sqrt(-2.*log(y));
00069     double p=t-((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1);
00070     return p;
00071   }
00072 }
00073 
00078 double inv_cumd_norm(const double& x)
00079 {
00080   double y=inv_cumd_norm_inner(x);
00081   y+=2.50662827*exp(.5*y*y)*(x-cumd_norm(y));
00082   return y;
00083 }
00084 
00090 double cumd_norm(const double& x)
00091 {
00092   const double b1=0.319381530;
00093   const double b2=-0.356563782;
00094   const double b3=1.781477937;
00095   const double b4=-1.821255978;
00096   const double b5=1.330274429;
00097   const double p=.2316419;
00098   if (x>=0)
00099   {
00100     double u=1./(1+p*x);
00101     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00102     double z=1.0-0.3989422804*exp(-.5*x*x)*y;
00103     return z;
00104   }
00105   else
00106   {
00107     double w=-x;
00108     double u=1./(1+p*w);
00109     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00110     double z=0.3989422804*exp(-.5*x*x)*y;
00111     return z;
00112   }
00113 }
00114 
00119 double bounded_cumd_norm(const double x,double beta)
00120 {
00121   const double b1=0.319381530;
00122   const double b2=-0.356563782;
00123   const double b3=1.781477937;
00124   const double b4=-1.821255978;
00125   const double b5=1.330274429;
00126   const double p=.2316419;
00127   if (x>=0)
00128   {
00129     double u=1./(1+p*x);
00130     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00131     double z=1.0-0.3989422804*exp(-.5*x*x)*y;
00132     z=beta*(z-0.5)+0.5;
00133     return z;
00134   }
00135   else
00136   {
00137     double w=-x;
00138     double u=1./(1+p*w);
00139     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00140     double z=0.3989422804*exp(-.5*x*x)*y;
00141     z=beta*(z-0.5)+0.5;
00142     return z;
00143   }
00144 }