ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
normmix.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: normmix.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 static double cc=0.39894228040143267794;   // 1/sqrt(2*pi)
00014 
00015 typedef double (*pinit_f)(double y,double a);
00016 
00021 static double cumd_normal_mixture(double x,double a)
00022 {
00023   // "normal" value for a is 3.0
00024   double y=0.95*cumd_norm(x)+0.05*cumd_norm(x/a);
00025   return y;
00026 }
00027 
00032 static double df_cumd_normal_mixture(double x,double a)
00033 {
00034   double x2=x*x;
00035   double dfx=cc*(0.95*exp(-0.5*x2)+0.05/a*exp(-0.5*x2/(a*a)) );
00036 
00037   return dfx;
00038 }
00039 
00044 static double cumd_normal_mixture_initx(double y,double a)
00045 {
00046   double x;
00047   if (y>0.999)
00048   {
00049     x= a*inv_cumd_norm((y-0.95)/0.05);
00050   }
00051   else if (y<.001)
00052   {
00053     x= 1.0-a*inv_cumd_norm((.05-y)/0.05);
00054   }
00055   else
00056   {
00057     x=inv_cumd_norm(y);
00058   }
00059   return x;
00060 }
00061 
00066 double  nr_generic(double y,pinit_f p_get_initial_x,
00067   pinit_f pfun,pinit_f pdfun)
00068 {
00069   double x=(*p_get_initial_x)(y,3.0);
00070 
00071   const int imax=15;
00072   int icount=0;
00073   do
00074   {
00075     icount++;
00076     double cy=(*pfun)(x,3.0);
00077     double der=(*pdfun)(x,3.0);
00078     double diff=y-cy;
00079     double h=diff/der;
00080     x+=h;
00081     if (fabs(h)<1.e-12) break;
00082   }
00083   while(icount<imax);
00084   //cout << " x = " << x << " icount = " << icount << endl;
00085 
00086 
00087   return x;
00088 }
00089 
00094 dvariable inv_cumd_normal_mixture(const prevariable& _yy,double a)
00095 {
00096   ADUNCONST(dvariable,yy)
00097   double  x=nr_generic(value(yy),cumd_normal_mixture_initx,cumd_normal_mixture,
00098     df_cumd_normal_mixture);
00099   dvariable z;
00100   value(z)=x;
00101   double dfx=1.0/df_cumd_normal_mixture(x,a);
00102   gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00103      &(value(z)), &(value(yy)),dfx);
00104   return z;
00105 }