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