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