ADMB Documentation  11.1.2403
 All Classes Files Functions Variables Typedefs Friends Defines
normmix2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: normmix2.cpp 1921 2014-04-23 18:47:05Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <fvar.hpp>
00012 
00013 typedef double (*pinit_f)(double y,double a);
00014 
00015 /*
00016 static double cc=0.39894228040143267794;
00017 static double cumd_normal_logistic_mixture(double x,double a)
00018 {
00019   // "normal" value for a is 3.0
00020   double y;
00021   if (x>-20.0)
00022   {
00023     y=0.95*cumd_norm(x)+0.05/(1.0+exp(-x/a));
00024   }
00025   else
00026   {
00027     y=0.95*cumd_norm(x)+0.05/(1.0+exp(-x/a));
00028   }
00029   return y;
00030 }
00031 
00032 static double df_cumd_normal_logistic_mixture(double x,double a)
00033 {
00034   // "normal" value for a is 3.0
00035   //double y=0.95*cumd_norm(x)+0.05*cumd_norm(x/a)
00036   double x2=x*x;
00037   double dfx;
00038   if (x>-20.0)
00039   {
00040     dfx=cc*0.95*exp(-0.5*x2)+0.05/a*exp(-x/a)/square(1.0+exp(-x/a));
00041   }
00042   else
00043   {
00044     dfx=cc*0.95*exp(-0.5*x2)+0.05/a*exp(-x/a)/square(1.0+exp(-x/a));
00045   }
00046 
00047   return dfx;
00048 }
00049 
00050 static double cumd_normal_logistic_mixture_initx(double y,double a)
00051 {
00052   double x;
00053   if (y>0.999)
00054   {
00055     x= a*inv_cumd_logistic((y-0.95)/0.05);
00056   }
00057   else if (y<.001)
00058   {
00059     x= 1.0-a*inv_cumd_logistic((.05-y)/0.05);
00060   }
00061   else
00062   {
00063     x=inv_cumd_norm(y);
00064   }
00065   return x;
00066 }
00067 */
00068 
00073 double  nr_generic(double y,double a,pinit_f p_get_initial_x,
00074   pinit_f pfun,pinit_f pdfun)
00075 {
00076   double x=(*p_get_initial_x)(y,a);
00077 
00078   const int imax=15;
00079   int icount=0;
00080   double h;
00081   do
00082   {
00083     icount++;
00084     double cy=(*pfun)(x,a);
00085     double der=(*pdfun)(x,a);
00086     double diff=y-cy;
00087     h=diff/der;
00088     x+=h;
00089     if (fabs(h)<1.e-12) break;
00090   }
00091   while(icount<imax);
00092   //cout << " x = " << x << " icount = " << icount << endl;
00093   if (fabs(h)>1.e-8)
00094   {
00095     cerr << "shit" << endl;
00096   }
00097 
00098   return x;
00099 }
00100 
00101 /*
00102 main()
00103 {
00104   gradient_structure gs(10000);
00105   dvariable y;
00106   cin >> y;
00107   dvariable x=inv_cumd_normal_logistic_mixture(y,3.0);
00108 }
00109 */