ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
cnorlogmix.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: cnorlogmix.cpp 1683 2014-02-26 00:16:59Z 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;
00014 
00015 typedef double (*pinit_f)(double y,double a);
00016 
00021 static double cumd_normal_logistic_mixture(double x,double a)
00022 {
00023   // "normal" value for a is 3.0
00024   double y;
00025   if (x>-20.0)
00026   {
00027     y=0.95*cumd_norm(x)+0.05/(1.0+exp(-x/a));
00028   }
00029   else
00030   {
00031     y=0.95*cumd_norm(x)+0.05/(1.0+exp(-x/a));
00032   }
00033   return y;
00034 }
00035 
00040 static double df_cumd_normal_logistic_mixture(double x,double a)
00041 {
00042   // "normal" value for a is 3.0
00043   //double y=0.95*cumd_norm(x)+0.05*cumd_norm(x/a)
00044   double x2=x*x;
00045   double dfx;
00046   if (x>-20.0)
00047   {
00048     dfx=cc*0.95*exp(-0.5*x2)+0.05/a*exp(-x/a)/square(1.0+exp(-x/a));
00049   }
00050   else
00051   {
00052     dfx=cc*0.95*exp(-0.5*x2)+0.05/a*exp(-x/a)/square(1.0+exp(-x/a));
00053   }
00054 
00055   return dfx;
00056 }
00057 
00062 static double cumd_normal_logistic_mixture_initx(double y,double a)
00063 {
00064   double x;
00065   if (y>0.999)
00066   {
00067     x= a*inv_cumd_logistic((y-0.95)/0.05);
00068   }
00069   else if (y<.001)
00070   {
00071     x= 1.0-a*inv_cumd_logistic((.05-y)/0.05);
00072   }
00073   else
00074   {
00075     x=inv_cumd_norm(y);
00076   }
00077   return x;
00078 }
00079 
00084 double  nr_generic(double y,double a,pinit_f p_get_initial_x,
00085   pinit_f pfun,pinit_f pdfun);
00086 
00087 dvariable inv_cumd_normal_logistic_mixture(const prevariable& _yy,double a)
00088 {
00089   ADUNCONST(dvariable,yy)
00090   double  x=nr_generic(value(yy),a,cumd_normal_logistic_mixture_initx,
00091     cumd_normal_logistic_mixture,df_cumd_normal_logistic_mixture);
00092   dvariable z;
00093   value(z)=x;
00094   double dfx=1.0/df_cumd_normal_logistic_mixture(x,a);
00095   gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00096      &(value(z)), &(value(yy)),dfx);
00097   return z;
00098 }
00099 
00100 /*
00101 main()
00102 {
00103   gradient_structure gs(10000);
00104   dvariable y;
00105   cin >> y;
00106   dvariable x=inv_cumd_normal_logistic_mixture(y,3.0);
00107 }
00108 */
00109 //********************************************************
00110 //********************************************************
00111 //********************************************************
00112 #include <fvar.hpp>
00113 
00118 double robust_normal_logistic_mixture_deviate(double x,double spread)
00119 {
00120   double y=cumd_norm(x);
00121   y = 0.99999999*y + 0.000000005; // To gain numerical stability
00122   double z = inv_cumd_normal_logistic_mixture(y,spread);
00123   return z;
00124 }
00125 
00130 double robust_normal_mixture_deviate(double x,double spread)
00131 {
00132   double y=cumd_norm(x);
00133   y = 0.99999999*y + 0.000000005; // To gain numerical stability
00134   double z = inv_cumd_normal_mixture(y,spread);
00135   return z;
00136 }
00137 
00142 double inv_cumd_normal_logistic_mixture(double yy,double a)
00143 {
00144   double  x=nr_generic(yy,a,cumd_normal_logistic_mixture_initx,
00145     cumd_normal_logistic_mixture,df_cumd_normal_logistic_mixture);
00146   return x;
00147 }