ADMB Documentation  11.1.2428
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2normmix2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2normmix2.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 
00013 static double cc=0.39894228040143267794;
00014 
00015 typedef double (*pinit_f)(double y,double a);
00016 
00017 double  nr_generic(double y,double a,pinit_f p_get_initial_x,
00018   pinit_f pfun,pinit_f pdfun);
00019 
00024 static double cumd_normal_logistic_mixture(double x,double a)
00025 {
00026   // "normal" value for a is 3.0
00027   double y;
00028   y=0.95*cumd_norm(x)+0.05/(1.0+exp(-x/a));
00029  /*
00030   if (x>-20.0)
00031   {
00032     y=0.95*cumd_norm(x)+0.05/(1.0+exp(-x));
00033   }
00034   else
00035   {
00036     y=0.95*cumd_norm(x)+0.05*exp(x)/(1.0+exp(x));
00037   }
00038   */
00039   return y;
00040 }
00041 
00046 static double df_cumd_normal_logistic_mixture(double x,double a)
00047 {
00048   // "normal" value for a is 3.0
00049   //double y=0.95*cumd_norm(x)+0.05*cumd_norm(x/a)
00050   double x2=x*x;
00051   double dfx;
00052     dfx=cc*0.95*exp(-0.5*x2)+0.05/a*exp(-x/a)/square(1.0+exp(-x/a));
00053  /*
00054   if (x>-20.0)
00055   {
00056     dfx=cc*0.95*exp(-0.5*x2)+0.05*exp(-x)/square(1.0+exp(-x));
00057   }
00058   else
00059   {
00060     dfx=cc*0.95*exp(-0.5*x2)+0.05*exp(x)/square(1.0+exp(x));
00061   }
00062   */
00063 
00064   return dfx;
00065 }
00066 
00071 static double cumd_normal_logistic_mixture_initx(double y,double a)
00072 {
00073   double x;
00074   if (y>0.999)
00075   {
00076     x= a*inv_cumd_logistic((y-0.95)/0.05);
00077   }
00078   else if (y<.001)
00079   {
00080     x= 1.0-a*inv_cumd_logistic((.05-y)/0.05);
00081   }
00082   else
00083   {
00084     x=inv_cumd_norm(y);
00085   }
00086   return x;
00087 }
00088 
00093 df1b2variable inv_cumd_normal_logistic_mixture(const df1b2variable& _yy,
00094   double a)
00095 {
00096   ADUNCONST(df1b2variable,yy)
00097   df1b2variable z;
00098   double  x=nr_generic(value(yy),a,cumd_normal_logistic_mixture_initx,
00099     cumd_normal_logistic_mixture,
00100     df_cumd_normal_logistic_mixture);
00101 
00102   // the derivatives of the inverse
00103   double x2=x*x;
00104   double e1=exp(-0.5*x2);
00105 
00106   double dgx=cc*0.95*e1+0.05/a*exp(-x/a)/square(1.0+exp(-x/a));
00107 
00108  /*
00109   double xp=x+1.e-5;
00110   double xm=x-1.e-5;
00111   double e1p=exp(-0.5*xp*xp);
00112   double e1m=exp(-0.5*xm*xm);
00113   double dgxp=cc*0.95*e1p+0.05*exp(-xp)/square(1.0+exp(-xp));
00114   double dgxm=cc*0.95*e1m+0.05*exp(-xm)/square(1.0+exp(-xm));
00115  */
00116 
00117   //cout << (dgxp-dgxm)/(2.0*1.e-5) << endl;
00118   //cout << (dgxp-2.0*dgx+dgxm)/(1.e-10) << endl;
00119 
00120   double d2g=-cc*0.95*x*e1+0.05/(a*a)*(
00121            -exp(-x/a)/square(1.0+exp(-x/a))
00122            +2.0*exp(-2.0*x/a)/cube(1.0+exp(-x/a)));
00123 
00124   double d3g=-cc*0.95*e1 +cc*x2*0.95*e1 +0.05/(a*a*a)*(
00125      exp(-x/a)/square(1.0+exp(-x/a)) -6.0*exp(-2.0*x/a)/cube(1.0+exp(-x/a))
00126      +6.0*exp(-3.0*x/a)/square(square(1.0+exp(-x/a))));
00127 
00128   double dfx=1.0/dgx;
00129   double d2f=-d2g*cube(dfx);
00130   double d3f=-d3g*cube(dfx)*dfx-3.0*d2g*d2f*square(dfx);
00131 
00132   double * yd=yy.get_u_dot();
00133   double * zd=z.get_u_dot();
00134   *z.get_u()=x;
00135   for (int i=0;i<df1b2variable::nvar;i++)
00136   {
00137     *zd++ =dfx * *yd++;
00138   }
00139   if (!df1b2_gradlist::no_derivatives)
00140     f1b2gradlist->write_pass1(&yy,&z,dfx,d2f,d3f);
00141 
00142   return z;
00143 }