ADMB Documentation  11.1.1903
 All Classes Files Functions Variables Typedefs Friends Defines
gammdev.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: gammdev.cpp 1107 2013-07-11 21:56:15Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #define HOME_VERSION
00012 #include "df1b2fun.h"
00013   df3_two_variable cumd_exponential(const df3_two_variable& x,
00014     const df3_two_variable& a);
00015 
00016 df3_two_variable cumd_cauchy(const df3_two_variable& x,
00017   const df3_two_variable& a);
00018 
00019   df1b2variable inv_cumd_gamma(const df1b2variable& _y,const df1b2variable& _a);
00020   double inv_cumd_gamma(double y,double _a);
00021 
00022   df3_two_variable cumd_gamma(const df3_two_variable& x,
00023     const df3_two_variable& a);
00024 
00029   df1b2variable gamma_deviate(const df1b2variable& _x,const df1b2variable& _a)
00030   {
00031     df1b2variable& x= (df1b2variable&)(_x);
00032     df1b2variable& a= (df1b2variable&)(_a);
00033 
00034     df1b2variable y=cumd_norm(x);
00035     y=.9999*y+.00005;
00036 
00037     //df1b2variable z=inv_cumd_gamma(y,a);
00038     df1b2variable z=inv_cumd_gamma(y,a);
00039 
00040     return z;
00041   }
00042 
00047   df1b2variable inv_cumd_gamma(const df1b2variable& _y,const df1b2variable& _a)
00048   {
00049     df1b2variable& y= (df1b2variable&)(_y);
00050     df1b2variable& a= (df1b2variable&)(_a);
00051     // get the inverse values
00052     double x=inv_cumd_gamma(value(y),value(_a));
00053 
00054     init_df3_two_variable xx(value(x));
00055     init_df3_two_variable aa(value(a));
00056     //init_df3_two_variable xx(1.0);
00057     //init_df3_two_variable aa(2.0);
00058     *xx.get_u_x()=1.0;
00059     *aa.get_u_y()=1.0;
00060 
00061     df3_two_variable z=cumd_gamma(xx,aa);
00062 
00063     // now use the derivatives of z to get the
00064     //derivatives of x wrt y,a and save them
00065 
00066     //double ca=value(a);
00067 
00068     double F_x=1.0/(*z.get_u_x());
00069 
00070     double F_y=-F_x* *z.get_u_y();
00071     double F_xx=-F_x* *z.get_u_xx()/square(*z.get_u_x());
00072 
00073     double F_xy=-(F_xx * *z.get_u_x() * *z.get_u_y() + F_x * *z.get_u_xy())/
00074                 *z.get_u_x();
00075 
00076     double F_yy=-(F_xx * square(*z.get_u_y())
00077                 +2.0* F_xy* *z.get_u_y()
00078                 +F_x * *z.get_u_yy());
00079 
00080     double F_xxx=-(F_x* *z.get_u_xxx()
00081                 +3.0*F_xx* *z.get_u_x() * *z.get_u_xx())
00082                 /cube(*z.get_u_x());
00083 
00084     double F_xxy=-(F_xxx * square(*z.get_u_x())* *z.get_u_y()
00085                  + 2.0*F_xx* *z.get_u_x()* *z.get_u_xy()
00086                  + F_xx* *z.get_u_xx() * *z.get_u_y()
00087                  + F_xy * *z.get_u_xx()
00088                  + F_x * *z.get_u_xxy())/
00089                  square(*z.get_u_x());
00090 
00091     double F_xyy=-(F_xxx* *z.get_u_x() *square(*z.get_u_y())
00092                  +2.0* F_xx * *z.get_u_xy() * *z.get_u_y()
00093                  +2.0*F_xxy * *z.get_u_x() * *z.get_u_y()
00094                  + 2.0*F_xy * *z.get_u_xy()
00095                  + F_xx * *z.get_u_x() * *z.get_u_yy()
00096                  + F_x * *z.get_u_xyy())/
00097                  *z.get_u_x();
00098      double F_yyy=-(F_xxx*cube(*z.get_u_y())+
00099                  +3.0*F_xxy*square(*z.get_u_y())
00100                  +3.0*F_xx* *z.get_u_y() * *z.get_u_yy()
00101                  +3.0*F_xy* *z.get_u_yy()
00102                  +3.0*F_xyy * *z.get_u_y()
00103                  +F_x * *z.get_u_yyy());
00104 
00105      df1b2variable zz;
00106      double * xd=_y.get_u_dot();
00107      double * yd=_a.get_u_dot();
00108      double * zd=zz.get_u_dot();
00109      *zz.get_u()=x;
00110 
00111      for (int i=0;i<df1b2variable::nvar;i++)
00112      {
00113        *zd++ = F_x * *xd++ + F_y * *yd++;
00114      }
00115      if (!df1b2_gradlist::no_derivatives)
00116      {
00117        f1b2gradlist->write_pass1(&_y,&_a,&zz,
00118         F_x,F_y,
00119         F_xx,F_xy,F_yy,
00120         F_xxx,F_xxy,F_xyy,F_yyy);
00121      }
00122      return zz;
00123   }