ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
prmonte.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: prmonte.cpp 2732 2014-11-26 22:54:51Z 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 dvector new_probing_bounded_multivariate_normal(int nvar, const dvector& a1,
00015   const dvector& b1, dmatrix& ch, const double& _wght,double pprobe,
00016   random_number_generator& rng)
00017 {
00018   double& wght=(double&) _wght;
00019   const double rob1=0.95;
00020   const double sqrt_tpi =sqrt(2*PI);
00021   dvector w(1,nvar);
00022   dvector a(1,nvar);
00023   dvector b(1,nvar);
00024   dvector alpha(1,nvar);
00025   dvector beta(1,nvar);
00026   a=a1;
00027   b=b1;
00028   wght=0;
00029   double rwght=0;
00030   double nwght=0;
00031   w.initialize();
00032   double ah;
00033   double bl;
00034   double upper;
00035   double lower;
00036   double upper1;
00037   double lower1;
00038   double diff;
00039   double diff1;
00040   int expflag;
00041   double y;
00042   //int in=0;
00043   //int ie=0;
00044   double _u = rng.better_rand();
00045   int rflag;
00046   if (_u>rob1)
00047   {
00048     rflag=1;
00049   }
00050   else
00051   {
00052     rflag=0;
00053   }
00054   if (!rflag)
00055   {
00056     for (int i=1;i<=nvar;i++)
00057     {
00058       ah=a(i)/ch(i,i);
00059       bl=b(i)/ch(i,i);
00060       double u = rng.better_rand();
00061       upper=cumd_norm(bl);
00062       lower=cumd_norm(ah);
00063       diff=upper-lower;
00064       if (diff>1.e-5)
00065       {
00066         expflag=0;
00067       }
00068       else
00069       {
00070         expflag=1;
00071       }
00072       upper1=cumd_cauchy(bl);
00073       lower1=cumd_cauchy(ah);
00074       diff1=upper1-lower1;
00075 
00076       u=u*.9998+.0001;
00077       if (!expflag)
00078       {
00079         y = inv_cumd_norm(u*(upper-lower)+lower);
00080       }
00081       else
00082       {
00083         y = inv_cumd_cauchy(u*(upper1-lower1)+lower1);
00084       }
00085       if (diff>1.e-5)
00086       {
00087         nwght+=-.5*y*y -log(sqrt_tpi*diff);
00088       }
00089       else
00090       {
00091         nwght += log_density_cauchy(y)-log(diff1);
00092       }
00093       for (int j=i;j<=nvar;j++)
00094       {
00095         double tmp=y*ch(j,i);
00096         w(j)+=tmp;
00097         a(j)-=tmp;
00098         b(j)-=tmp;
00099       }
00100     }
00101     wght = nwght;
00102   }
00103   else
00104   {
00105     a=a1;
00106     b=b1;
00107     for (int i=1;i<=nvar;i++)
00108     {
00109       ah=a(i)/ch(i,i);
00110       bl=b(i)/ch(i,i);
00111       double u = rng.better_rand();
00112       double pp = rng.better_rand();
00113       upper=cumd_norm(bl);
00114       lower=cumd_norm(ah);
00115       diff=upper-lower;
00116       if (diff>1.e-5)
00117       {
00118         expflag=0;
00119       }
00120       else
00121       {
00122         expflag=1;
00123       }
00124       upper1=cumd_cauchy(bl);
00125       lower1=cumd_cauchy(ah);
00126       diff1=upper1-lower1;
00127 
00128       u=u*.9998+.0001;
00129       if (!expflag)
00130       {
00131         if (pp>pprobe)
00132           y = inv_cumd_norm(u*(upper-lower)+lower);
00133         else
00134           y = inv_cumd_cauchy(u*(upper1-lower1)+lower1);
00135       }
00136       else
00137       {
00138         y = inv_cumd_cauchy(u*(upper1-lower1)+lower1);
00139       }
00140       if (diff>1.e-5)
00141       {
00142         rwght+=log((1.0-pprobe)*exp(-.5*y*y)/(sqrt_tpi*diff)
00143          +pprobe*density_cauchy(y)/diff1);
00144       }
00145       else
00146       {
00147         rwght += log_density_cauchy(y)-log(diff1);
00148       }
00149       for (int j=i;j<=nvar;j++)
00150       {
00151         double tmp=y*ch(j,i);
00152         w(j)+=tmp;
00153         a(j)-=tmp;
00154         b(j)-=tmp;
00155       }
00156     }
00157     double dd=rob1*(exp(nwght-rwght))+1.0-rob1;
00158     if (dd<=0)
00159     {
00160       cerr << "dd <=0" << endl;
00161     }
00162     wght = log(dd)+rwght;
00163   }
00164   return w;
00165 }
00166 
00167 void new_probing_bounded_multivariate_normal_mcmc(int nvar, const dvector& a1,
00168   const dvector& b1, dmatrix& ch, const double& _wght, const dvector& _y,
00169   double pprobe, random_number_generator& rng)
00170 {
00171   double& wght=(double&) _wght;
00172   const double rob1=0.95;
00173   const double sqrt_tpi =sqrt(2*PI);
00174   dvector w(1,nvar);
00175   dvector a(1,nvar);
00176   dvector b(1,nvar);
00177   dvector alpha(1,nvar);
00178   dvector beta(1,nvar);
00179   wght=0;
00180   double rwght=0;
00181   double nwght=0;
00182   w.initialize();
00183   double ah;
00184   double bl;
00185   double upper;
00186   double lower;
00187   double upper1;
00188   double lower1;
00189   double diff;
00190   double diff1;
00191   int expflag;
00192   double y = 0;
00193   //int in=0;
00194   //int ie=0;
00195 
00196   a=a1;
00197   b=b1;
00198   int i;
00199   for (i=1;i<=nvar;i++)
00200   {
00201     y=_y(i);
00202     ah=a(i)/ch(i,i);
00203     bl=b(i)/ch(i,i);
00204     /*double u = */rng.better_rand();
00205     upper=cumd_norm(bl);
00206     lower=cumd_norm(ah);
00207     diff=upper-lower;
00208 /*
00209     if (diff>1.e-5)
00210     {
00211       expflag=0;
00212     }
00213     else
00214     {
00215       expflag=1;
00216     }
00217 */
00218     upper1=cumd_cauchy(bl);
00219     lower1=cumd_cauchy(ah);
00220     diff1=upper1-lower1;
00221 
00222     if (diff>1.e-5)
00223     {
00224       nwght+=-.5*y*y -log(sqrt_tpi*diff);
00225     }
00226     else
00227     {
00228       nwght+= log_density_cauchy(y)-log(diff1);
00229     }
00230     for (int j=i;j<=nvar;j++)
00231     {
00232       double tmp=y*ch(j,i);
00233       w(j)+=tmp;
00234       a(j)-=tmp;
00235       b(j)-=tmp;
00236     }
00237   }
00238   a=a1;
00239   b=b1;
00240   w.initialize();
00241   for (i=1;i<=nvar;i++)
00242   {
00243     //y=_y(i);
00244     ah=a(i)/ch(i,i);
00245     bl=b(i)/ch(i,i);
00246     double u = rng.better_rand();
00247     double pp = rng.better_rand();
00248     upper=cumd_norm(bl);
00249     lower=cumd_norm(ah);
00250     diff=upper-lower;
00251     if (diff>1.e-5)
00252     {
00253       expflag=0;
00254     }
00255     else
00256     {
00257       expflag=1;
00258     }
00259     upper1=cumd_cauchy(bl);
00260     lower1=cumd_cauchy(ah);
00261     diff1=upper1-lower1;
00262 
00263     u=u*.9998+.0001;
00264     if (!expflag)
00265     {
00266       if (pp>pprobe)
00267         y = inv_cumd_norm(u*(upper-lower)+lower);
00268       else
00269         y = inv_cumd_cauchy(u*(upper1-lower1)+lower1);
00270     }
00271     else
00272     {
00273       y = inv_cumd_cauchy(u*(upper1-lower1)+lower1);
00274     }
00275     if (diff>1.e-5)
00276     {
00277       rwght+=log((1.0-pprobe)*exp(-.5*y*y)/(sqrt_tpi*diff)
00278        +pprobe*density_cauchy(y)/diff1);
00279     }
00280     else
00281     {
00282       rwght += log_density_cauchy(y)-log(diff1);
00283     }
00284   }
00285   wght = log(rob1*(exp(nwght-rwght))+(1.0-rob1))+rwght;
00286   for (int j=i;j<=nvar;j++)
00287   {
00288     double tmp=y*ch(j,i);
00289     w(j)+=tmp;
00290     a(j)-=tmp;
00291     b(j)-=tmp;
00292   }
00293 }