ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2normmix.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2normmix.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,pinit_f p_get_initial_x,
00018   pinit_f pfun,pinit_f pdfun);
00019 
00024 static double cumd_normal_mixture(double x,double a)
00025 {
00026   double y=0.95*cumd_norm(x)+0.05*cumd_norm(x/a);
00027   return y;
00028 }
00029 
00034 static double df_cumd_normal_mixture(double x,double a)
00035 {
00036   double x2=x*x;
00037   double dfx=cc*(0.95*exp(-0.5*x2)+0.05/a*exp(-0.5*x2/(a*a)) );
00038 
00039   return dfx;
00040 }
00041 
00046 static double cumd_normal_mixture_initx(double y,double a)
00047 {
00048   double x;
00049   if (y>0.999)
00050   {
00051     x= a*inv_cumd_norm((y-0.95)/0.05);
00052   }
00053   else if (y<.001)
00054   {
00055     x= 1.0-a*inv_cumd_norm((.05-y)/0.05);
00056   }
00057   else
00058   {
00059     x=inv_cumd_norm(y);
00060   }
00061   return x;
00062 }
00063 
00068 df1b2variable inv_cumd_normal_mixture(const df1b2variable& _yy,double a)
00069 {
00070   ADUNCONST(df1b2variable,yy)
00071   df1b2variable z;
00072   double  x=nr_generic(value(yy),cumd_normal_mixture_initx,cumd_normal_mixture,
00073     df_cumd_normal_mixture);
00074 
00075   // the derivatives of the inverse
00076   double x2=x*x;
00077   double a2=a*a;
00078   double a3=a2*a;
00079   double a5=a2*a3;
00080   double e1=exp(-0.5*x2);
00081   double e2=exp(-0.5*x2/a2);
00082 
00083   double dgx=cc*(0.95*e1+0.05/a*e2);
00084 
00085   double d2g=-cc*x*(0.95*e1+0.05/a3*e2);
00086 
00087   double d3g=-cc*(0.95*e1+0.05/a3*e2)
00088     +cc*x2*(0.95*e1+0.05/a5*e2);
00089 
00090   double dfx=1.0/dgx;
00091   double d2f=-d2g*cube(dfx);
00092   double d3f=-d3g*cube(dfx)*dfx-3.0*d2g*d2f*square(dfx);
00093 
00094   double * yd=yy.get_u_dot();
00095   double * zd=z.get_u_dot();
00096   *z.get_u()=x;
00097   for (int i=0;i<df1b2variable::nvar;i++)
00098   {
00099     *zd++ =dfx * *yd++;
00100   }
00101   if (!df1b2_gradlist::no_derivatives)
00102     f1b2gradlist->write_pass1(&yy,&z,dfx,d2f,d3f);
00103 
00104   return z;
00105 }
00106