ADMB Documentation  11.1.1913
 All Classes Files Functions Variables Typedefs Friends Defines
mod_mc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_mc.cpp 1827 2014-04-03 20:24:43Z 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 inv_cumd_cauchy(const double& x);
00011 double inv_cumd_norms(const double& x);
00012 double cumd_norm(const double& x);
00013 double cumd_cauchy(const double& x);
00014 double density_cauchy(const double& x);
00015 double myran1(long int&);
00016 double better_rand(long int&);
00017 
00018 double log_likelihood_mixture(const double& x);
00019 
00020 void multivariate_mixture(const dvector& _mix,int nvar,long int& iseed,
00021   const double& _log_density_normal, const double& _log_density_cauchy,
00022   const double& _log_density_small_normal,int is);
00023 
00024 dvector bounded_multivariate_cauchy(int nvar, const dvector& a1,
00025   dvector& b1, const dmatrix& ch,long int& iseed, const double& lprob,
00026   double& log_tprob, const int& outflag);
00027 
00028 dvector bounded_robust_multivariate_normal(int nvar, const dvector& a1,
00029   dvector& b1, const dmatrix& ch, const dmatrix& ch3, const dmatrix& chinv,
00030   const dmatrix& ch3inv, double contaminant,long int& iseed,
00031   const double& lprob, const double& lprob3, double& log_tprob,
00032   const int& outflag);
00033 
00034 void function_minimizer::monte_carlo_routine(void)
00035 {
00036   initial_params::mc_phase=1;
00037   if (stddev_params::num_stddev_params==0) return;
00038   {
00039     int nvar=initial_params::nvarcalc(); // get the number of active parameters
00040     dvector x(1,nvar);
00041     dvector jac(1,nvar);
00042     initial_params::xinit(x);
00043     initial_params::stddev_scale(jac,x);
00044     dvector bmn(1,nvar);
00045     bmn.initialize();
00046 // ***************************************************************
00047     dvector scale(1,nvar);   // need to get scale from somewhere
00048     dvector diag(1,nvar);   // need to get scale from somewhere
00049     dvector v(1,nvar);  // need to read in v from model.rep
00050     dmatrix S(1,nvar,1,nvar);
00051     {
00052       adstring tmpstring="admodel.cov";
00053       if (ad_comm::wd_flag)
00054         tmpstring = ad_comm::adprogram_name + ".cov";
00055       uistream cif((char*)tmpstring);
00056       if (!cif)
00057       {
00058         cerr << "Error trying to open file " << tmpstring
00059              << " for reading" << endl;
00060       }
00061       int tmp_nvar = 0;
00062       cif >> tmp_nvar;
00063       if (nvar !=tmp_nvar)
00064       {
00065         cerr << "Incorrect number of independent variables in file"
00066             << tmpstring  << endl;
00067         exit(1);
00068       }
00069       cif >> S;
00070       if (!cif)
00071       {
00072         cerr << "error reading covariance matrix from "
00073              <<  tmpstring << endl;
00074         exit(1);
00075       }
00076 
00077       initial_params::mc_phase=0;
00078       /*int check=*/initial_params::stddev_scale(scale,x);
00079       initial_params::mc_phase=1;
00080       {
00081         dmatrix tmp(1,nvar,1,nvar);
00082 
00083         int i;
00084         for (i=1;i<=nvar;i++)
00085         {
00086           tmp(i,i)=S(i,i)*scale(i)*scale(i);
00087           for (int j=1;j<i;j++)
00088           {
00089             tmp(i,j)=S(i,j)*scale(i)*scale(j);
00090             tmp(j,i)=tmp(i,j);
00091           }
00092           diag(i)=sqrt(tmp(i,i));
00093         }
00094         S=tmp;
00095         //cout << endl << S << endl << endl;
00096         //cout << endl << choleski_decomp(S) << endl << endl;
00097 
00098         for (i=1;i<=nvar;i++)
00099         {
00100           S(i,i)=S(i,i)/(diag(i)*diag(i));
00101           for (int j=1;j<i;j++)
00102           {
00103             S(i,j)=S(i,j)/(diag(i)*diag(j));
00104             S(j,i)=S(i,j);
00105           }
00106         }
00107       }
00108     }
00109     dmatrix CHD= choleski_decomp(S);
00110     int sgn;
00111     /*double lnd=*/ln_det(S,sgn);
00112     //cout << endl << S << endl << endl;
00113     //cout << endl << CHD << endl << endl;
00114 
00115     dmatrix chdec(1,nvar,1,nvar);
00116     chdec=choleski_decomp(S);
00117     {
00118       ofstream ofs("chol2");
00119       ofs << chdec << endl;
00120     }
00121     //cout << endl << chdec << endl << endl;
00122 
00123     //int check=initial_params::stddev_scale(jac,x);
00124     initial_params::montecarlo_scale(scale,x);
00125     {
00126       ofstream ofs("admodel.tmp");
00127       dmatrix covar(1,nvar,1,nvar);
00128     }
00129     double ll=0.0;
00130     initial_params::add_random_vector(jac,bmn,ll,diag);
00131 /*
00132 #if defined(USE_DDOUBLE)
00133     double ljac0=sum(log(jac+double(1.e-50)));
00134 #else
00135     double ljac0=sum(log(jac+1.e-50));
00136 #endif
00137 */
00138 
00139     dmatrix ch3=3.*chdec;
00140     dmatrix chinv=inv(chdec);
00141     dmatrix ch3inv=inv(ch3);
00142     //double xxin;
00143 
00144     {
00145       long int iseed=0;
00146       int number_sims=  500;
00147       //cin >> iseed;
00148       iseed=-33517;
00149       cout << "Enter value for seed" << endl;
00150       cin >> iseed;
00151       cout << "Enter number of simulations" << endl;
00152       cin >> number_sims;
00153       if (iseed>0)
00154       {
00155         iseed=-iseed;
00156       }
00157       better_rand(iseed);
00158       //double lprob=0.0;
00159       //double lprob3=0.0;
00160       // get lower and upper bounds
00161 
00162       //independent_variables x(1,nvar);
00163       independent_variables parsave(1,nvar);
00164 
00165       int ii=1;
00166       initial_params::copy_all_values(parsave,ii);
00167       gradient_structure::set_NO_DERIVATIVES();
00168       //ofstream ooff((char*) (ad_comm::adprogram_name + adstring(".mte")) );
00169       //ofstream ooff1((char*) (ad_comm::adprogram_name + adstring(".mt2")) );
00170       //double cont=0.00;
00171       double log_tprob_normal=0.0;
00172       double log_tprob_small_normal=0.0;
00173       double log_tprob_cauchy=0.0;
00174       //double log_tprob=0.0;
00175       int ndvar=stddev_params::num_stddev_calc();
00176       dvector param_values(1,ndvar);
00177       //int outflag;
00178       //ooff1 << "Number of simulations = " << number_sims << endl;
00179       ii=1;
00180       initial_params::restore_all_values(parsave,ii);
00181 
00182       //double fbest;
00183 
00184       if (!ad_comm::pvm_manager)
00185       {
00186         /*fbest=*/get_monte_carlo_value(nvar,x);
00187       }
00188       else
00189       {
00190         switch (ad_comm::pvm_manager->mode)
00191         {
00192         case 1: // master
00193           /*fbest=*/pvm_master_get_monte_carlo_value(nvar,x);
00194           break;
00195         case 2: // slave
00196           pvm_slave_get_monte_carlo_value(nvar);
00197           break;
00198         default:
00199           cerr << "error illegal value for pvm_manager->mode" << endl;
00200           exit(1);
00201         }
00202       }
00203 
00204       multivariate_mixture(bmn,nvar,iseed,log_tprob_normal,
00205         log_tprob_cauchy,log_tprob_small_normal,-1);
00206       bmn=elem_div(bmn,scale);
00207 /*
00208       if (log_tprob_normal >= log_tprob_cauchy)
00209       {
00210         log_tprob=log_tprob_normal
00211           +log(0.95+0.05*exp(log_tprob_cauchy-log_tprob_normal));
00212       }
00213       else
00214       {
00215         log_tprob=log_tprob_cauchy
00216           +log(0.95*exp(log_tprob_normal-log_tprob_cauchy)+.05);
00217       }
00218 */
00219       //double cdiff=-fbest-log_tprob;
00220       //double cfb=-fbest;
00221       //double clt=log_tprob;
00222       //ooff1 << " *  weight    likelihood   simprob  ln det "
00223        //"  ljac0   ljac     parameter value " << endl;
00224       //ooff1 << setw(10) << exp(-fbest-log_tprob) << " "
00225         //<< setw(10) << exp(-fbest) << " "
00226       //  << setw(10) << exp(log_tprob) << " "
00227        //    << setw(10) << lnd << " "
00228        //    << setw(10) << param_values << endl;
00229 
00230       dvector y(1,nvar);
00231 
00232       // what is this supposed to do?
00233       //initial_params::get_jacobian_value(y,jac);
00234       //for (int i=1;i<=nvar;i++)
00235       //{
00236        // for (int j=1;j<=nvar;j++)
00237         //{
00238          // chdec(i,j)/=jac(i);
00239         //}
00240       //}
00241 
00242       //ofstream ogs("sims");
00243       //ogs << nvar << " " << number_sims << endl;
00244       for (int i=1;i<=number_sims;i++)
00245       {
00246         log_tprob_normal=0.0;
00247         log_tprob_small_normal=0.0;
00248         log_tprob_cauchy=0.0;
00249         //log_tprob=0.0;
00250         ii=1;
00251         initial_params::restore_all_values(parsave,ii);
00252 
00253         double mixprob=better_rand(iseed);
00254         int mixswitch;
00255         //if (mixprob<.0)
00256         if (mixprob<.05)
00257         {
00258           mixswitch=1;
00259         }
00260         else if (mixprob<.50)
00261         {
00262           mixswitch=2;
00263         }
00264         else
00265         {
00266           mixswitch=0;
00267         }
00268         multivariate_mixture(bmn,nvar,iseed,log_tprob_normal,
00269           log_tprob_cauchy,log_tprob_small_normal,mixswitch);
00270         //bmn=elem_div(bmn,scale);
00271 
00272 /*
00273         if (log_tprob_normal >= log_tprob_cauchy)
00274         {
00275           log_tprob=log_tprob_normal
00276             +log( 0.50+.45*exp(log_tprob_small_normal-log_tprob_normal)
00277             +      .05*exp(log_tprob_cauchy-log_tprob_normal));
00278         }
00279         else
00280         {
00281           log_tprob=log_tprob_cauchy
00282             +log( 0.50*exp(log_tprob_normal-log_tprob_cauchy)
00283             +     0.45*exp(log_tprob_small_normal-log_tprob_cauchy)+.05 );
00284         }
00285 */
00286         dvector bmn1=CHD*bmn;
00287         //bmn1=elem_div(bmn1,scale);
00288 
00289         ll=0.0;
00290         initial_params::add_random_vector(jac,bmn1,ll,diag);
00291         initial_params::xinit(y);
00292 
00293         //initial_params::stddev_scale(jac,y);
00294 /*
00295 #if defined(USE_DDOUBLE)
00296         double ljac=sum(log(jac+double(1.e-50)));
00297 #else
00298         double ljac=sum(log(jac+1.e-50));
00299 #endif
00300 */
00301 
00302 
00303        // ogs << log_tprob << " " << ll << " " << x << endl;
00304         if (!ad_comm::pvm_manager)
00305         {
00306           get_monte_carlo_value(nvar,x);
00307         }
00308         else
00309         {
00310           switch (ad_comm::pvm_manager->mode)
00311           {
00312           case 1: // master
00313             pvm_master_get_monte_carlo_value(nvar,x);
00314             break;
00315           case 2: // slave
00316             pvm_master_get_monte_carlo_value(nvar,x);
00317             break;
00318           default:
00319             cerr << "error illega value for pvm_manager->mode" << endl;
00320             exit(1);
00321           }
00322         }
00323 
00324         //ooff << setw(12) << -f-log_tprob << "  "
00325         //     << setw(12) << -f-log_tprob-ll << "  "
00326         //     << setw(12) << fbest-f << "  "
00327         //     << setw(12) << log_tprob << " "
00328         //     << setw(5) << outflag << " " << setw(12) << -f;
00329         //ooff << endl;
00330 
00331         ii=1;
00332         stddev_params::copy_all_values(param_values,ii);
00333         //ooff << " " << setw(6) << "  " << param_values << endl;
00334         /*
00335         if (mixswitch==1)
00336         {
00337           ooff1 << " *  weight    likelihood   simprob  ln det "
00338              "  ljac0   ljac     parameter value " << endl;
00339         }
00340         else if (mixswitch==2)
00341         {
00342           ooff1 << "    weight    likelihood   simprob  ln det "
00343              "  ljac0   ljac     parameter value " << endl;
00344         }
00345         else
00346         {
00347           ooff1 << " +  weight    likelihood   simprob  ln det "
00348              "  ljac0   ljac     parameter value " << endl;
00349         }
00350         ooff1 << setw(10) << exp(-f-log_tprob-cdiff+(ljac0-ljac) ) << " "
00351           << setw(10) << -f << " "
00352           << setw(10) << log_tprob << " "
00353           << setw(10) << lnd << " "
00354           << setw(10) << ljac0 << " "
00355           << setw(10) << ljac << " "
00356           << setw(10) << param_values << endl;
00357         */
00358       }
00359       initial_params::mc_phase=0;
00360     }
00361   }
00362 }