ADMB Documentation  11.1.2551
 All Classes Files Functions Variables Typedefs Friends Defines
prmonte.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: prmonte.cpp 1964 2014-04-30 22:09:22Z 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     if (diff>1.e-5)
00209     {
00210       expflag=0;
00211     }
00212     else
00213     {
00214       expflag=1;
00215     }
00216     upper1=cumd_cauchy(bl);
00217     lower1=cumd_cauchy(ah);
00218     diff1=upper1-lower1;
00219 
00220     if (diff>1.e-5)
00221     {
00222       nwght+=-.5*y*y -log(sqrt_tpi*diff);
00223     }
00224     else
00225     {
00226       nwght+= log_density_cauchy(y)-log(diff1);
00227     }
00228     for (int j=i;j<=nvar;j++)
00229     {
00230       double tmp=y*ch(j,i);
00231       w(j)+=tmp;
00232       a(j)-=tmp;
00233       b(j)-=tmp;
00234     }
00235   }
00236   a=a1;
00237   b=b1;
00238   w.initialize();
00239   for (i=1;i<=nvar;i++)
00240   {
00241     y=_y(i);
00242     ah=a(i)/ch(i,i);
00243     bl=b(i)/ch(i,i);
00244     double u = rng.better_rand();
00245     double pp = rng.better_rand();
00246     upper=cumd_norm(bl);
00247     lower=cumd_norm(ah);
00248     diff=upper-lower;
00249     if (diff>1.e-5)
00250     {
00251       expflag=0;
00252     }
00253     else
00254     {
00255       expflag=1;
00256     }
00257     upper1=cumd_cauchy(bl);
00258     lower1=cumd_cauchy(ah);
00259     diff1=upper1-lower1;
00260 
00261     u=u*.9998+.0001;
00262     if (!expflag)
00263     {
00264       if (pp>pprobe)
00265         y = inv_cumd_norm(u*(upper-lower)+lower);
00266       else
00267         y = inv_cumd_cauchy(u*(upper1-lower1)+lower1);
00268     }
00269     else
00270     {
00271       y = inv_cumd_cauchy(u*(upper1-lower1)+lower1);
00272     }
00273     if (diff>1.e-5)
00274     {
00275       rwght+=log((1.0-pprobe)*exp(-.5*y*y)/(sqrt_tpi*diff)
00276        +pprobe*density_cauchy(y)/diff1);
00277     }
00278     else
00279     {
00280       rwght += log_density_cauchy(y)-log(diff1);
00281     }
00282   }
00283   wght = log(rob1*(exp(nwght-rwght))+(1.0-rob1))+rwght;
00284   for (int j=i;j<=nvar;j++)
00285   {
00286     double tmp=y*ch(j,i);
00287     w(j)+=tmp;
00288     a(j)-=tmp;
00289     b(j)-=tmp;
00290   }
00291 }