ADMB Documentation  11.1.2497
 All Classes Files Functions Variables Typedefs Friends Defines
xmonte2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xmonte2.cpp 1919 2014-04-22 22:02:01Z 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 inv_cumd_norm(const double& x);
00010 double cumd_norm(const double& x);
00011 double myran1(long int&);
00012 //double better_rand(long int&);
00013 
00014 void bounded_multivariate_normal_mcmc(int nvar, const dvector& a1,
00015   const dvector& b1, dmatrix& ch, const double& _wght, const dvector& y,
00016   const random_number_generator& rng)
00017 {
00018   double & wght=(double &) _wght;
00019   //cout << y << endl;
00020   const double sqrt_tpi =sqrt(2*PI);
00021   dvector a(1,nvar);
00022   dvector b(1,nvar);
00023   dvector alpha(1,nvar);
00024   dvector beta(1,nvar);
00025   a=a1;
00026   b=b1;
00027   wght=0;
00028   double ah;
00029   double bl;
00030   double upper;
00031   double lower;
00032   double diff;
00033   int expflag;
00034   int in=0;
00035   int ie=0;
00036   for (int i=1;i<=nvar;i++)
00037   {
00038     ah=a(i)/ch(i,i);
00039     bl=b(i)/ch(i,i);
00040     upper=cumd_norm(bl);
00041     lower=cumd_norm(ah);
00042     diff=upper-lower;
00043     if (diff>1.e-5)
00044     {
00045       wght-=log(diff);
00046       expflag=0;
00047     }
00048     else
00049     {
00050       upper=cumd_cauchy(bl);
00051       lower=cumd_cauchy(ah);
00052       diff=upper-lower;
00053       wght-=log(diff);
00054       expflag=1;
00055     }
00056 
00057     if (!expflag)
00058     {
00059       wght -= .5*y(i)*y(i);
00060       in++;
00061     }
00062     else
00063     {
00064       ie++;
00065       wght += log_density_cauchy(y(i));
00066     }
00067 
00068 
00069     for (int j=i;j<=nvar;j++)
00070     {
00071       double tmp=y(i)*ch(j,i);
00072       a(j)-=tmp;
00073       b(j)-=tmp;
00074     }
00075   }
00076   wght +=  in*log(1./sqrt_tpi);
00077 }
00078 
00079 void probing_bounded_multivariate_normal_mcmc(int nvar, const dvector& a1,
00080   const dvector& b1, dmatrix& ch, const double& _wght, const dvector& y,
00081   double pprobe, const random_number_generator& rng)
00082 {
00083   double & wght=(double &) _wght;
00084   //cout << y << endl;
00085   const double sqrt_tpi =sqrt(2*PI);
00086   dvector a(1,nvar);
00087   dvector b(1,nvar);
00088   dvector alpha(1,nvar);
00089   dvector beta(1,nvar);
00090   a=a1;
00091   b=b1;
00092   wght=0;
00093   double ah;
00094   double bl;
00095   double upper;
00096   double lower;
00097   double diff;
00098   double diff1;
00099   //int in=0;
00100   //int ie=0;
00101   for (int i=1;i<=nvar;i++)
00102   {
00103     ah=a(i)/ch(i,i);
00104     bl=b(i)/ch(i,i);
00105     upper=cumd_norm(bl);
00106     lower=cumd_norm(ah);
00107     diff=upper-lower;
00108     upper=cumd_cauchy(bl);
00109     lower=cumd_cauchy(ah);
00110     diff1=upper-lower;
00111     if (diff>1.e-5)
00112     {
00113       wght+=log((1.0-pprobe)*exp(-.5*y(i)*y(i))/(sqrt_tpi*diff)
00114          +pprobe*density_cauchy(y(i))/diff1);
00115     }
00116     else
00117     {
00118       wght += log_density_cauchy(y(i))-log(diff1);
00119     }
00120 
00121     for (int j=i;j<=nvar;j++)
00122     {
00123       double tmp=y(i)*ch(j,i);
00124       a(j)-=tmp;
00125       b(j)-=tmp;
00126     }
00127   }
00128 }
00129 
00130 void bounded_multivariate_uniform_mcmc(int nvar, const dvector& a1,
00131   const dvector& b1, dmatrix& ch, const double& _wght, const dvector& y,
00132   const random_number_generator& rng)
00133 {
00134   double& wght=(double&) _wght;
00135   dvector a(1,nvar);
00136   dvector b(1,nvar);
00137   a=a1;
00138   b=b1;
00139   wght=0;
00140   double ah;
00141   double bl;
00142   double upper;
00143   double lower;
00144   double diff;
00145   for (int i=1;i<=nvar;i++)
00146   {
00147     ah=a(i)/ch(i,i);
00148     bl=b(i)/ch(i,i);
00149     lower=ffmax(-1.0,ah);
00150     upper=ffmin(1.0,bl);
00151     diff=upper-lower;
00152     wght-=log(diff);
00153     for (int j=i;j<=nvar;j++)
00154     {
00155       double tmp=y(i)*ch(j,i);
00156       a(j)-=tmp;
00157       b(j)-=tmp;
00158     }
00159   }
00160 }