ADMB Documentation  11.1x.2735
 All Classes Files Functions Variables Typedefs Friends Defines
mod_mc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_mc.cpp 1971 2014-05-02 18:45:14Z 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 #ifdef USE_ADPVM
00185       if (ad_comm::pvm_manager)
00186       {
00187         switch (ad_comm::pvm_manager->mode)
00188         {
00189         case 1: // master
00190           /*fbest=*/pvm_master_get_monte_carlo_value(nvar,x);
00191           break;
00192         case 2: // slave
00193           pvm_slave_get_monte_carlo_value(nvar);
00194           break;
00195         default:
00196           cerr << "error illegal value for pvm_manager->mode" << endl;
00197           exit(1);
00198         }
00199       }
00200       else
00201 #endif
00202       {
00203         /*fbest=*/get_monte_carlo_value(nvar,x);
00204       }
00205 
00206       multivariate_mixture(bmn,nvar,iseed,log_tprob_normal,
00207         log_tprob_cauchy,log_tprob_small_normal,-1);
00208       bmn=elem_div(bmn,scale);
00209 /*
00210       if (log_tprob_normal >= log_tprob_cauchy)
00211       {
00212         log_tprob=log_tprob_normal
00213           +log(0.95+0.05*exp(log_tprob_cauchy-log_tprob_normal));
00214       }
00215       else
00216       {
00217         log_tprob=log_tprob_cauchy
00218           +log(0.95*exp(log_tprob_normal-log_tprob_cauchy)+.05);
00219       }
00220 */
00221       //double cdiff=-fbest-log_tprob;
00222       //double cfb=-fbest;
00223       //double clt=log_tprob;
00224       //ooff1 << " *  weight    likelihood   simprob  ln det "
00225        //"  ljac0   ljac     parameter value " << endl;
00226       //ooff1 << setw(10) << exp(-fbest-log_tprob) << " "
00227         //<< setw(10) << exp(-fbest) << " "
00228       //  << setw(10) << exp(log_tprob) << " "
00229        //    << setw(10) << lnd << " "
00230        //    << setw(10) << param_values << endl;
00231 
00232       dvector y(1,nvar);
00233 
00234       // what is this supposed to do?
00235       //initial_params::get_jacobian_value(y,jac);
00236       //for (int i=1;i<=nvar;i++)
00237       //{
00238        // for (int j=1;j<=nvar;j++)
00239         //{
00240          // chdec(i,j)/=jac(i);
00241         //}
00242       //}
00243 
00244       //ofstream ogs("sims");
00245       //ogs << nvar << " " << number_sims << endl;
00246       for (int i=1;i<=number_sims;i++)
00247       {
00248         log_tprob_normal=0.0;
00249         log_tprob_small_normal=0.0;
00250         log_tprob_cauchy=0.0;
00251         //log_tprob=0.0;
00252         ii=1;
00253         initial_params::restore_all_values(parsave,ii);
00254 
00255         double mixprob=better_rand(iseed);
00256         int mixswitch;
00257         //if (mixprob<.0)
00258         if (mixprob<.05)
00259         {
00260           mixswitch=1;
00261         }
00262         else if (mixprob<.50)
00263         {
00264           mixswitch=2;
00265         }
00266         else
00267         {
00268           mixswitch=0;
00269         }
00270         multivariate_mixture(bmn,nvar,iseed,log_tprob_normal,
00271           log_tprob_cauchy,log_tprob_small_normal,mixswitch);
00272         //bmn=elem_div(bmn,scale);
00273 
00274 /*
00275         if (log_tprob_normal >= log_tprob_cauchy)
00276         {
00277           log_tprob=log_tprob_normal
00278             +log( 0.50+.45*exp(log_tprob_small_normal-log_tprob_normal)
00279             +      .05*exp(log_tprob_cauchy-log_tprob_normal));
00280         }
00281         else
00282         {
00283           log_tprob=log_tprob_cauchy
00284             +log( 0.50*exp(log_tprob_normal-log_tprob_cauchy)
00285             +     0.45*exp(log_tprob_small_normal-log_tprob_cauchy)+.05 );
00286         }
00287 */
00288         dvector bmn1=CHD*bmn;
00289         //bmn1=elem_div(bmn1,scale);
00290 
00291         ll=0.0;
00292         initial_params::add_random_vector(jac,bmn1,ll,diag);
00293         initial_params::xinit(y);
00294 
00295         //initial_params::stddev_scale(jac,y);
00296 /*
00297 #if defined(USE_DDOUBLE)
00298         double ljac=sum(log(jac+double(1.e-50)));
00299 #else
00300         double ljac=sum(log(jac+1.e-50));
00301 #endif
00302 */
00303 
00304 
00305        // ogs << log_tprob << " " << ll << " " << x << endl;
00306 #ifdef USE_ADPVM
00307         if (ad_comm::pvm_manager)
00308         {
00309           switch (ad_comm::pvm_manager->mode)
00310           {
00311           case 1: // master
00312             pvm_master_get_monte_carlo_value(nvar,x);
00313             break;
00314           case 2: // slave
00315             pvm_master_get_monte_carlo_value(nvar,x);
00316             break;
00317           default:
00318             cerr << "error illega value for pvm_manager->mode" << endl;
00319             exit(1);
00320           }
00321         }
00322         else
00323 #endif
00324         {
00325           get_monte_carlo_value(nvar,x);
00326         }
00327 
00328         //ooff << setw(12) << -f-log_tprob << "  "
00329         //     << setw(12) << -f-log_tprob-ll << "  "
00330         //     << setw(12) << fbest-f << "  "
00331         //     << setw(12) << log_tprob << " "
00332         //     << setw(5) << outflag << " " << setw(12) << -f;
00333         //ooff << endl;
00334 
00335         ii=1;
00336         stddev_params::copy_all_values(param_values,ii);
00337         //ooff << " " << setw(6) << "  " << param_values << endl;
00338         /*
00339         if (mixswitch==1)
00340         {
00341           ooff1 << " *  weight    likelihood   simprob  ln det "
00342              "  ljac0   ljac     parameter value " << endl;
00343         }
00344         else if (mixswitch==2)
00345         {
00346           ooff1 << "    weight    likelihood   simprob  ln det "
00347              "  ljac0   ljac     parameter value " << endl;
00348         }
00349         else
00350         {
00351           ooff1 << " +  weight    likelihood   simprob  ln det "
00352              "  ljac0   ljac     parameter value " << endl;
00353         }
00354         ooff1 << setw(10) << exp(-f-log_tprob-cdiff+(ljac0-ljac) ) << " "
00355           << setw(10) << -f << " "
00356           << setw(10) << log_tprob << " "
00357           << setw(10) << lnd << " "
00358           << setw(10) << ljac0 << " "
00359           << setw(10) << ljac << " "
00360           << setw(10) << param_values << endl;
00361         */
00362       }
00363       initial_params::mc_phase=0;
00364     }
00365   }
00366 }