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