ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
mod_mc3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_mc3.cpp 1942 2014-04-28 22:22:45Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 double better_rand(long int&);
00010 
00011 void initial_params::add_random_vector(const dvector& y, const dvector& x,
00012  const double& ll, const dvector& diag)
00013 {
00014   int ii=1;
00015   for (int i=0;i<num_initial_params;i++)
00016   {
00017     if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00018     {
00019       (varsptr[i])->add_value(y,x,ii,ll,diag);
00020     }
00021   }
00022 }
00023 void initial_params::get_jacobian_value(const dvector& y, const dvector& jac)
00024 {
00025   int ii=1;
00026   for (int i=0;i<num_initial_params;i++)
00027   {
00028     if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00029     {
00030       (varsptr[i])->get_jacobian(y,jac,ii);
00031     }
00032   }
00033 }
00034 
00035 void multivariate_mixture(const dvector& _mix, int nvar, long int& iseed,
00036   const double& _log_density_normal, const double& _log_density_cauchy,
00037   const double& _log_density_small_normal, int is)
00038 {
00039   dvector& mix=(dvector&) _mix;
00040   double& log_density_cauchy=(double&) _log_density_cauchy;
00041   double& log_density_normal=(double&) _log_density_normal;
00042   double& log_density_small_normal=(double&) _log_density_small_normal;
00043   const double r2=sqrt(2.0);
00044   const double l2p=0.5*log(2*PI);
00045   const double l3p=0.5*log(2*PI)-log(3.0);
00046   const double pr2=log(PI*r2);
00047   log_density_normal=0.0;
00048   log_density_cauchy=0.0;
00049 
00050   if (is==0)
00051   {
00052     for (int i=1;i<=nvar;i++)
00053     {
00054       double u;
00055       do
00056       {
00057         u = better_rand(iseed);
00058       }
00059       while (u<.0001 || u>.9999);
00060       mix(i) = inv_cumd_norm(u);
00061       log_density_normal-= l2p +.5*mix(i)*mix(i);
00062       log_density_small_normal-= l3p +4.5*mix(i)*mix(i);
00063       log_density_cauchy+= -pr2 + 1./(1.+mix(i)*mix(i));
00064     }
00065   }
00066   else if (is==2)
00067   {
00068     for (int i=1;i<=nvar;i++)
00069     {
00070       double u;
00071       do
00072       {
00073         u = better_rand(iseed);
00074       }
00075       while (u<.0001 || u>.9999);
00076       mix(i) = inv_cumd_norm(u);
00077       mix(i)/=3;
00078       log_density_normal-= l2p +.5*mix(i)*mix(i);
00079       log_density_cauchy+= -pr2 + 1./(1.+mix(i)*mix(i));
00080       log_density_small_normal-= l3p +4.5*mix(i)*mix(i);
00081     }
00082   }
00083   else if (is==1)
00084   {
00085     for (int i=1;i<=nvar;i++)
00086     {
00087       double u;
00088       do
00089       {
00090         u = better_rand(iseed);
00091       }
00092       while (u<.0001 || u>.9999);
00093       mix(i) = inv_cumd_cauchy(u);
00094       log_density_normal-= l2p +.5*mix(i)*mix(i);
00095       log_density_small_normal-= l3p +4.5*mix(i)*mix(i);
00096       log_density_cauchy+= -pr2 + 1./(1.+mix(i)*mix(i));
00097     }
00098   }
00099   else
00100   {
00101     for (int i=1;i<=nvar;i++)
00102     {
00103       double u=0.5;
00104       while (u<.0001 || u>.9999);
00105       mix(i) = inv_cumd_cauchy(u);
00106       log_density_normal-= l2p +.5*mix(i)*mix(i);
00107       log_density_small_normal-= l3p +4.5*mix(i)*mix(i);
00108       log_density_cauchy+= -pr2 + 1./(1.+mix(i)*mix(i));
00109     }
00110   }
00111 }
00112 
00113 double set_value_inv_mc(const prevariable& z,double min, double max)
00114 {
00115   return tan(PI*( (value(z)-min)/(max-min)-0.5 ));
00116 }
00117 
00118 double set_value_inv_mc(double z,double min, double max)
00119 {
00120   return tan(PI*( (z-min)/(max-min)-0.5 ));
00121 }
00122 
00123 dvariable set_value_mc(const prevariable& z,double min,double max)
00124 {
00125   const double pinv=1./PI;
00126   dvariable y=atan(z)*pinv+0.5;
00127   return min+(max-min)*y;
00128 }
00129 
00130 double set_value_mc(double z,double min,double max)
00131 {
00132   const double pinv=1./PI;
00133   double y=atan(z)*pinv+0.5;
00134   return min+(max-min)*y;
00135 }
00136 
00137 void set_value_inv_mc(const dvar_vector& x, const dvector& _v, const int& _ii,
00138   const double fmin, const double fmax)
00139 {
00140   dvector& v=(dvector&) _v;
00141   int& ii=(int&) _ii;
00142   int min=x.indexmin();
00143   int max=x.indexmax();
00144   for (int i=min;i<=max;i++)
00145   {
00146     v(ii++)=set_value_inv_mc(x(i),fmin,fmax);
00147   }
00148 }
00149 
00150 void set_value_mc(const dvar_vector& _x,const dvar_vector& v, const int& _ii,
00151   const double fmin, const double fmax)
00152 {
00153   ADUNCONST(dvar_vector,x)
00154   int& ii=(int&) _ii;
00155   int min=x.indexmin();
00156   int max=x.indexmax();
00157   for (int i=min;i<=max;i++)
00158   {
00159     x(i)=set_value_mc(v(ii++),fmin,fmax);
00160   }
00161 }