ADMB Documentation: vcumdist.cpp (* $Id$) Source File

Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <fvar.hpp>
00012 
00013 dvariable inv_cumd_norm(const prevariable& x);
00014 prevariable& cumd_norm(const prevariable& x);
00015 /*
00016 double normal_tail_right(const double& x)
00017 {
00018   const double a3=5;
00019   const double a4=9;
00020   const double a5=129;
00021   double x2=x*x;
00022   double z=0.3989422804*exp(-.5*x2);
00023   double b1=x2+2;
00024   double b2=b1*(x2+4);
00025   double b3=b2*(x2+6);
00026   double b4=b3*(x2+8);
00027   double b5=b4*(x2+10);
00028   double tmp=z/x*(1.0 -1.0/b1 +1.0/b2 - a3/b3 +a4/b4 -a5/b5);
00029   return tmp;
00030 }
00031 */
00032 
00037 dvariable inv_cumd_norm_inner(const prevariable& x)
00038 {
00039  if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN)
00040    gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
00041  RETURN_ARRAYS_INCREMENT();
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_ARRAYS_DECREMENT();
00052     return 0.0;
00053   }
00054 
00055   if (x<=0.5)
00056   {
00057     //dvariable t = sqrt(-2.*log(x));
00058     //dvariable p=((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1)-t;
00059     //RETURN_ARRAYS_DECREMENT();
00060     //return p;
00061 
00062     double tt = sqrt(-2.*log(value(x)));
00063     double u=(c2*tt+c1)*tt+c0;
00064     double v=((d3*tt+d2)*tt+d1)*tt+1;
00065     double vinv=1.0/v;
00066     double pp=u*vinv-tt;
00067 
00068     //double pp=u*vinv-tt;
00069     double dfu=vinv;
00070     double dfvinv=u;
00071     double dftt=-1.0;
00072     //double vinv=1.0/v;
00073     double dfv=-vinv*vinv*dfvinv;
00074     //double v=((d3*tt+d2)*tt+d1)*tt+1;
00075     dftt+=((3*d3*tt+2.0*d2)*tt+d1)*dfv;
00076     //double u=(c2*tt+c1)*tt+c0;
00077     dftt+=(2.0*c2*tt+c1)*dfu;
00078     //double tt = sqrt(-2.*log(value(x)));
00079     double dfx=-1.0/(tt*value(x))*dftt;
00080 
00081     RETURN_ARRAYS_DECREMENT();
00082     gradient_structure::RETURN_PTR->v->x=pp;
00083     gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00084        &(gradient_structure::RETURN_PTR->v->x), &(x.v->x),dfx);
00085     return(*gradient_structure::RETURN_PTR);
00086   }
00087   else if (x==0.5)
00088   {
00089     cout << "can't happen" << endl;
00090     exit(1);
00091     //return 0.0;
00092   }
00093   else
00094   {
00095     //dvariable y=1.-x;
00096     //dvariable t = sqrt(-2.*log(y));
00097 
00098     //dvariable p=t-((c2*t+c1)*t+c0)/((((d3*t+d2)*t+d1)*t)+1);
00099     //RETURN_ARRAYS_DECREMENT();
00100     //return p;
00101 
00102     double yy=1.-value(x);
00103     double tt = sqrt(-2.*log(yy));
00104     double u=((c2*tt+c1)*tt+c0);
00105     double v=((d3*tt+d2)*tt+d1)*tt+1;
00106     double vinv=1/v;
00107     double pp=tt-u*vinv;
00108 
00109     //double pp=tt-u*vinv;
00110     double dfu=-vinv;
00111     double dfvinv=-u;
00112     double dftt=1.0;
00113     //double vinv=1.0/v;
00114     double dfv=-vinv*vinv*dfvinv;
00115     //double v=((d3*tt+d2)*tt+d1)*tt+1;
00116     dftt+=((3*d3*tt+2.0*d2)*tt+d1)*dfv;
00117     //double u=(c2*tt+c1)*tt+c0;
00118     dftt+=(2.0*c2*tt+c1)*dfu;
00119     //double tt = sqrt(-2.*log(yy));
00120     double dfy=-1.0/(tt*yy)*dftt;
00121     //double yy=1.-value(x);
00122     double dfx=-dfy;
00123 
00124     RETURN_ARRAYS_DECREMENT();
00125     gradient_structure::RETURN_PTR->v->x=pp;
00126     gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00127        &(gradient_structure::RETURN_PTR->v->x), &(x.v->x),dfx);
00128     return(*gradient_structure::RETURN_PTR);
00129   }
00130 }
00131 
00136 dvariable inv_cumd_norm(const prevariable& x)
00137 {
00138   dvariable y=inv_cumd_norm_inner(x);
00139   if (x>1.e-30)
00140     y+=2.50662827*exp(.5*y*y)*(x-cumd_norm(y));
00141   return y;
00142 }
00143 
00148 dvariable old_cumd_norm(const prevariable& x)
00149 {
00150   RETURN_ARRAYS_INCREMENT();
00151   const double b1=0.319381530;
00152   const double b2=-0.356563782;
00153   const double b3=1.781477937;
00154   const double b4=-1.821255978;
00155   const double b5=1.330274429;
00156   const double p=.2316419;
00157   if (x>=0)
00158   {
00159     dvariable u=1./(1+p*x);
00160     dvariable y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00161     dvariable z=1.0-0.3989422804*exp(-.5*x*x)*y;
00162     RETURN_ARRAYS_DECREMENT();
00163     return z;
00164   }
00165   else
00166   {
00167     dvariable w=-x;
00168     dvariable u=1./(1+p*w);
00169     dvariable y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00170     dvariable z=0.3989422804*exp(-.5*x*x)*y;
00171     RETURN_ARRAYS_DECREMENT();
00172     return z;
00173   }
00174 }
00175 
00181 prevariable& cumd_norm(const prevariable& _x)
00182 {
00183  if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN)
00184    gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
00185 
00186   double x=value(_x);
00187   const double b1=0.319381530;
00188   const double b2=-0.356563782;
00189   const double b3=1.781477937;
00190   const double b4=-1.821255978;
00191   const double b5=1.330274429;
00192   const double b55=b5*5;
00193   const double b44=b4*4;
00194   const double b33=b3*3;
00195   const double b22=b2*2;
00196   const double p=.2316419;
00197   if (x>=0)
00198   {
00199     double u=1./(1+p*x);
00200     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00201     double tmp1=-0.3989422804*exp(-.5*x*x);
00202     double z=1.0+tmp1*y;
00203     gradient_structure::RETURN_PTR->v->x=z;
00204 
00205 
00206     // double z=1.0+tmp1*y;
00207     double dftmp1=y;
00208     double dfy=tmp1;
00209 
00210     // double tmp1=-0.3989422804*exp(-.5*x*x);
00211     double dfx=-tmp1*x*dftmp1;
00212 
00213     //double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00214     double dfu=  ((((b55*u+b44)*u+b33)*u+b22)*u+b1)*dfy;
00215 
00216     //double u=1./(1+p*x);
00217     dfx-=u*u*p*dfu;
00218 
00219     gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00220        &(gradient_structure::RETURN_PTR->v->x), &(_x.v->x), dfx);
00221   }
00222   else
00223   {
00224     double w=-x;
00225     double u=1./(1+p*w);
00226     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00227     double tmp1=0.3989422804*exp(-.5*x*x);
00228     double z=tmp1*y;
00229 
00230     //double z=tmp1*y;
00231     double dftmp1=y;
00232     double dfy=tmp1;
00233 
00234     //double tmp1=0.3989422804*exp(-.5*x*x);
00235     double dfx=-tmp1*x*dftmp1;
00236 
00237     //double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00238     double dfu=  ((((b55*u+b44)*u+b33)*u+b22)*u+b1)*dfy;
00239 
00240     //double u=1./(1+p*w);
00241     double dfw=-u*u*p*dfu;
00242 
00243     //double w=-value(x);
00244     dfx-=dfw;
00245 
00246     gradient_structure::RETURN_PTR->v->x=z;
00247     gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00248        &(gradient_structure::RETURN_PTR->v->x), &(_x.v->x),dfx);
00249   }
00250   return(*gradient_structure::RETURN_PTR);
00251 }
00252 
00257 dvar_vector inv_cumd_norm(const dvar_vector& x)
00258 {
00259   int mmin=x.indexmin();
00260   int mmax=x.indexmax();
00261   dvar_vector tmp(mmin,mmax);
00262   for (int i=mmin;i<=mmax;i++)
00263   {
00264     tmp(i)=inv_cumd_norm(x(i));
00265   }
00266   return tmp;
00267 }
00268 
00273 prevariable& bounded_cumd_norm(const prevariable& _x,double beta)
00274 {
00275  if (++gradient_structure::RETURN_PTR > gradient_structure::MAX_RETURN)
00276    gradient_structure::RETURN_PTR = gradient_structure::MIN_RETURN;
00277 
00278   double x=value(_x);
00279   const double b1=0.319381530;
00280   const double b2=-0.356563782;
00281   const double b3=1.781477937;
00282   const double b4=-1.821255978;
00283   const double b5=1.330274429;
00284   const double b55=b5*5;
00285   const double b44=b4*4;
00286   const double b33=b3*3;
00287   const double b22=b2*2;
00288   const double p=.2316419;
00289   if (x>=0)
00290   {
00291     double u=1./(1+p*x);
00292     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00293     double tmp1=-0.3989422804*exp(-.5*x*x);
00294     double z1=1.0+tmp1*y;
00295     double z=beta*(z1-0.5)+0.5;
00296     gradient_structure::RETURN_PTR->v->x=z;
00297 
00298     //double z=beta*(z1-0.5)+0.5
00299     double dfz1=beta;
00300 
00301     // double z1=1.0+tmp1*y;
00302     double dftmp1=y*dfz1;
00303     double dfy=tmp1*dfz1;
00304 
00305     // double tmp1=-0.3989422804*exp(-.5*x*x);
00306     double dfx=-tmp1*x*dftmp1;
00307 
00308     //double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00309     double dfu=  ((((b55*u+b44)*u+b33)*u+b22)*u+b1)*dfy;
00310 
00311     //double u=1./(1+p*x);
00312     dfx-=u*u*p*dfu;
00313 
00314     gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00315        &(gradient_structure::RETURN_PTR->v->x), &(_x.v->x), dfx);
00316   }
00317   else
00318   {
00319     double w=-x;
00320     double u=1./(1+p*w);
00321     double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00322     double tmp1=0.3989422804*exp(-.5*x*x);
00323     double z1=tmp1*y;
00324     double z=beta*(z1-0.5)+0.5;
00325 
00326     //double z=beta*(z1-0.5)+0.5;
00327     double dfz1=beta;
00328 
00329     //double z1=tmp1*y;
00330     double dftmp1=y*dfz1;
00331     double dfy=tmp1*dfz1;
00332 
00333     //double tmp1=0.3989422804*exp(-.5*x*x);
00334     double dfx=-tmp1*x*dftmp1;
00335 
00336     //double y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00337     double dfu=  ((((b55*u+b44)*u+b33)*u+b22)*u+b1)*dfy;
00338 
00339     //double u=1./(1+p*w);
00340     double dfw=-u*u*p*dfu;
00341 
00342     //double w=-value(x);
00343     dfx-=dfw;
00344 
00345     gradient_structure::RETURN_PTR->v->x=z;
00346     gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00347        &(gradient_structure::RETURN_PTR->v->x), &(_x.v->x),dfx);
00348   }
00349   return(*gradient_structure::RETURN_PTR);
00350 }
00351 
00356 dvariable inv_cumd_norm_logistic(const prevariable& x,double p)
00357 {
00358 #if !defined(OPT_LIB)
00359   if (0.0<p || 1.0>p)
00360   {
00361     cerr << "Error in dvariable inv_cumd_norm_logistic -- illegal p value = "
00362          << p << endl;
00363     exit(1);
00364   }
00365 #endif
00366   dvariable y=inv_cumd_norm_inner(x);
00367   y+=( (1.0-p)*2.50662827*exp(.5*y*y) +
00368     p*exp(-x)/square(1.0+exp(-x)) ) * (x-cumd_norm_logistic (y,p));
00369   return y;
00370 }
00371 
00376 prevariable& cumd_norm_logistic(const prevariable& _x,double p)
00377 {
00378   return (1.0-p)*cumd_norm(_x)+p*cumd_logistic(_x);
00379 }