ADMB Documentation  11.1.2494
 All Classes Files Functions Variables Typedefs Friends Defines
cnormmix.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: cnormmix.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 
00040 static double cumd_normal_mixture_initx(double y,double a)
00041 {
00042   double x;
00043   if (y>0.999)
00044   {
00045     x= a*inv_cumd_norm((y-0.95)/0.05);
00046   }
00047   else if (y<.001)
00048   {
00049     x= 1.0-a*inv_cumd_norm((.05-y)/0.05);
00050   }
00051   else
00052   {
00053     x=inv_cumd_norm(y);
00054   }
00055   return x;
00056 }
00057 
00062 static double  nr_generic(double y,pinit_f p_get_initial_x,
00063   pinit_f pfun,pinit_f pdfun)
00064 {
00065   double x=(*p_get_initial_x)(y,3.0);
00066 
00067   const int imax=15;
00068   int icount=0;
00069   do
00070   {
00071     icount++;
00072     double cy=(*pfun)(x,3.0);
00073     double der=(*pdfun)(x,3.0);
00074     double diff=y-cy;
00075     double h=diff/der;
00076     x+=h;
00077     if (fabs(h)<1.e-12) break;
00078   }
00079   while(icount<imax);
00080   //cout << " x = " << x << " icount = " << icount << endl;
00081 
00082 
00083   return x;
00084 }
00085 
00090 double inv_cumd_normal_mixture(double yy,double a)
00091 {
00092   double  x=nr_generic(yy,cumd_normal_mixture_initx,cumd_normal_mixture,
00093     df_cumd_normal_mixture);
00094   return x;
00095 }