ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
mod_pmin.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_pmin.cpp 2721 2014-11-25 19:46:50Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 void get_confidence_interval(const dvector& left_bd, const dvector& right_bd,
00010   dmatrix& ms, const dvector& xs, const dvector& siglevel,
00011   const int& level_index, dvector& xdist,int index);
00012 void get_onesided_intervals(const dvector& left_bd, const dvector& right_bd,
00013   dmatrix& ms, const dvector& xs, const dvector& siglevel,
00014   const int& level_index, dvector& xdist,int index);
00015 void report_confidence_limits(const ofstream& ofs3,int numsig_levels,
00016   dvector& siglevel, const dvector& left_bd, const dvector& right_bd);
00017 void report_onesided_confidence_limits(const ofstream& ofs3,int numsig_levels,
00018   dvector& siglevel, const dvector& left_bd, const dvector& right_bd,int ip);
00019 
00020 
00021 dmatrix trans(const dvector& x)
00022 {
00023   int mmin=x.indexmin();
00024   int mmax=x.indexmax();
00025   dmatrix tmp(mmin,mmax,1,1);
00026   for (int i=mmin;i<=mmax;i++)
00027   {
00028     tmp(i,1)=x(i);
00029   }
00030   return tmp;
00031 }
00032   double mult_factor(int j)
00033   {
00034     switch(j)
00035     {
00036       case 0:
00037         return .0;
00038       case 1:
00039         return .5;
00040       case 2:
00041         return 1.0;
00042       case 3:
00043         return 1.5;
00044       default:
00045         return 2.0;
00046     }
00047   }
00048 
00049   double trimax(double x,double y,double z);
00050 
00051 #if defined(USE_ADPVM)
00052   void function_minimizer::pvm_slave_likeprof_routine(void)
00053   {
00054     do
00055     {
00056       int prof_switch=get_int_from_master();
00057       if (!prof_switch) break;
00058       if (prof_switch !=3)
00059       {
00060         cerr << "Error in prof_switch " << prof_switch << endl;
00061         ad_exit(1);
00062       }
00063       int underflow_flag=get_int_from_master();
00064       pvm_slave_prof_minimize(underflow_flag);
00065     }
00066     while(1);
00067   }
00068 #endif
00069   void function_minimizer::likeprof_routine(double global_min)
00070   {
00071     int on1 = 0;
00072     int nopt = 0;
00073     if ( (on1=option_match(ad_comm::argc,ad_comm::argv,"-iprint",nopt))>-1)
00074     {
00075       if (!nopt)
00076       {
00077         cerr << "Usage -iprint option needs integer  -- ignored" << endl;
00078         iprint = 10;
00079       }
00080       else
00081       {
00082         int jj=atoi(ad_comm::argv[on1+1]);
00083         iprint = jj;
00084       }
00085     }
00086     else
00087     {
00088       iprint = 10;
00089     }
00090 
00091     dvector siglevel("{.90,.95,.975}");
00092     int num_pp=likeprof_params::likeprofptr[0]->get_stepnumber();
00093     {
00094       for (int ip=1;ip<likeprof_params::num_likeprof_params;ip++)
00095       {
00096         int sno=likeprof_params::likeprofptr[ip]->get_stepnumber();
00097         if (sno)
00098         {
00099           if (sno>num_pp) num_pp=sno;
00100         }
00101       }
00102     }
00103     initial_params::current_phase = initial_params::max_number_phases;
00104     // DF NOV 28 11
00105     if (random_effects_flag)
00106     {
00107      initial_params::set_inactive_only_random_effects();
00108     }
00109     int nvar=initial_params::nvarcalc();
00110     dvector xsave(1,nvar);
00111     {
00112       int ii=1;
00113       initial_params::copy_all_values(xsave,ii);
00114     }
00115     double old_value;  // this is where we were
00116     double new_value;  // this is where we want to go
00117     double fprof = 0.0;
00118     double current_value;
00119     dmatrix lprof(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00120     dmatrix ldet(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00121     dmatrix ln_det_proj_jac(0,likeprof_params::num_likeprof_params-1,
00122       -num_pp,num_pp);
00123     dmatrix pdet(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00124     dmatrix actual_value(0,likeprof_params::num_likeprof_params-1,
00125       -num_pp,num_pp);
00126     dvector all_values(-num_pp,num_pp);
00127     dvector all_num_sigs(-num_pp,num_pp);
00128     //dvector xxxtmp(-num_pp,num_pp);
00129     //d3_array hesses(-num_pp,num_pp,1,nvar,1,nvar);
00130     dmatrix lg_jacob(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00131     dmatrix lg_prjacob(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00132     //dmatrix xxtmp(-num_pp,num_pp,1,nvar);
00133     dvector xvector(1,nvar);
00134     dmatrix xmax(-num_pp,num_pp,1,nvar);  // this holds the conditional max
00135     dmatrix gprof(-num_pp,num_pp,1,nvar);  // this holds the conditional max
00136     dmatrix fgrads(-num_pp,num_pp,1,nvar);  // this holds the conditional max
00137     dmatrix xdist(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00138     int sign=0;
00139     double sigma;
00140     ofstream ofs2((char*) (ad_comm::adprogram_name + adstring(".prf")) );
00141     double udet =unrestricted_hess_determinant();
00142     const int offset=0;
00143     dvector xscale(1,nvar);   // need to get scale from somewhere
00144     dvector likeprof_save(0,likeprof_params::num_likeprof_params-1);
00145     dmatrix penalties(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00146     penalties.initialize();
00147 #ifdef CURVE_CORRECT
00148     int nlp=likeprof_params::num_likeprof_params;
00149     d3_array eigenvals(0,nlp-1,-num_pp,num_pp,1,nvar-1);
00150     d3_array curvcor(0,nlp-1,-num_pp,num_pp,1,nvar-1);
00151 #endif
00152 
00153     int ip;
00154     for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00155     {
00156       likeprof_save(ip)=likeprof_params::likeprofptr[ip]->get_value();
00157     }
00158     double final_weight;
00159     for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00160     {
00161       int sno=likeprof_params::likeprofptr[ip]->get_stepnumber();
00162       double snz=likeprof_params::likeprofptr[ip]->get_stepsize();
00163       if (sno)
00164       {
00165         num_pp=sno;
00166       }
00167       //snz (step size) must be greater than zero.
00168       //Check "void likeprof_params::set_stepsize(double x)"
00169       const double relsig = snz > 0 ? snz : 0.5;
00170 
00171       if (ip>0)
00172       {
00173         int ii=1;
00174         initial_params::restore_all_values(xsave,ii);
00175       }
00176       sigma=likeprof_params::likeprofptr[ip]->get_sigma(); // this is the
00177       if (sigma==0.0)
00178         cerr << "error standard dev of likeporf parameter is 0" << endl;
00179                                         // estimated sd
00180       old_value=likeprof_save(ip);
00181       old_value=old_value+offset*relsig*sigma;  // this is where we
00182       int bigint_flag=0;
00183       int bigint_flag1=0;
00184       double ldiff=0.0;
00185       double num_sigs;
00186       for (int i=1;i<=2;i++)  // go in positive and negative directions
00187       {
00188         num_sigs=0.0;
00189         bigint_flag=0;
00190         bigint_flag1=0;
00191         int underflow_flag=0;
00192         if (i>1) // get the parameter values at the global minimum
00193         {
00194           int ii=1;
00195           initial_params::restore_all_values(xsave,ii);
00196         }
00197         if (i==1)
00198         {
00199           sign=1;
00200         }
00201         else
00202         {
00203           sign=-1;
00204         }
00205         current_value=old_value;  // initialize at the minimum
00206         for (int j=0;j<=num_pp;j++)  // go in positive and negative directions
00207         {
00208           if (j!=0 || sign > 0)
00209           {
00210             if (bigint_flag==0)
00211             {
00212               num_sigs+=mult_factor(j)*relsig*sign;
00213               current_value+=mult_factor(j)*relsig*sign*sigma;
00214               new_value=current_value;
00215               // new_value=current_value+j*relsig*sign*sigma;
00216             }
00217             else
00218             {
00219               if (bigint_flag1==0)
00220               {
00221                 num_sigs+=1.5*relsig*sign;
00222                 current_value+=1.5*relsig*sign*sigma;
00223                 new_value=current_value;
00224               }
00225               else
00226               {
00227                 num_sigs+=2.5*relsig*sign;
00228                 current_value+=2.5*relsig*sign*sigma;
00229                 new_value=current_value;
00230               }
00231             }
00232 #if defined(USE_ADPVM)
00233           if (ad_comm::pvm_manager)
00234           {
00235             if (ad_comm::pvm_manager->mode==1) // master
00236             {
00237               send_int_to_slaves(3);
00238               pvm_master_prof_minimize(ip,sigma,new_value,fprof,underflow_flag,
00239                 global_min,penalties(ip,j*sign),final_weight); // get the
00240             }
00241           }
00242           else
00243 #endif
00244           {
00245             if (random_effects_flag==0)
00246             {
00247               prof_minimize(ip,sigma,new_value,fprof,underflow_flag,
00248               global_min,penalties(ip,j*sign),final_weight); // get the
00249                                                         // conditional max
00250             }
00251             else
00252             {
00253               prof_minimize_re(ip,sigma,new_value,fprof,underflow_flag,
00254               global_min,penalties(ip,j*sign),final_weight); // get the
00255                                                         // conditional max
00256             }
00257           }
00258           all_num_sigs(j*sign)=num_sigs;
00259           initial_params::xinit(xvector); // save the
00260           int ic=1;
00261           // save the conditional maximum
00262           initial_params::copy_all_values(xmax(sign*j),ic);
00263           /*int check=*/initial_params::stddev_scale(xscale,xvector);
00264         //#if defined(DO_PROFILE)
00265           //dvector curvscale(1,nvar);   // need to get scale from somewhere
00266         //#endif
00267 
00268        // #if defined(DO_PROFILE)
00269           // check=initial_params::stddev_curvscale(curvscale,xvector);
00270        // #endif
00271           //cout << "xscale = " << endl << xscale << endl;
00272           //{
00273            // ofstream ofs("xscale");
00274            // ofs << xscale << endl;
00275          // }
00276 #if defined(USE_DDOUBLE)
00277           lg_jacob(ip,sign*j)=sum(log(fabs(xscale)+double(1.e-60)));
00278 #else
00279           lg_jacob(ip,sign*j)=sum(log(fabs(xscale)+1.e-60));
00280 #endif
00281           if (!underflow_flag)
00282           {
00283             lprof(ip,sign*j)=fprof;
00284           }
00285           else
00286           {
00287             lprof(ip,sign*j)=lprof(ip,sign*(j-1))+2;
00288           }
00289           double xx=likeprof_params::likeprofptr[ip]->get_value();
00290           if (!underflow_flag)
00291           {
00292             actual_value(ip,sign*j)=xx;
00293           }
00294           else
00295           {
00296             actual_value(ip,sign*j)=new_value;  // this is where we
00297           }
00298 
00299           ldiff=fprof-lprof(ip,0);
00300           if ( ldiff > 40.0)
00301           {
00302             underflow_flag=1;
00303           }
00304 
00305           //if ( (fprof-lprof(ip,0)) > 3.0 && bigint_flag ==0)
00306           if (abs(j) > 4)
00307           {
00308             bigint_flag=1;
00309           }
00310           if (abs(j) > 5)
00311           {
00312             bigint_flag1=1;
00313           }
00314           // get the gradient for the profile likelihood variable at
00315           // the conditional maximum
00316           initial_params::current_phase = initial_params::max_number_phases;
00317           nvar=initial_params::nvarcalc();
00318           dvector g(1,nvar);
00319           dvector fg(1,nvar);
00320           if (!underflow_flag)
00321           {
00322             // g is the grad. wrt the prof lik var
00323             // fg is the grad. wrt the log-likelihood
00324             get_particular_grad(ip,nvar,fg,g);
00325             gprof(sign*j)=g;
00326 
00327             fgrads(sign*j)=fg;
00328             xdist(ip,sign*j)=norm(elem_div(gprof(sign*j),xscale));
00329           }
00330           else
00331           {
00332             gprof(sign*j)=gprof(sign*(j-1));
00333             xdist(ip,sign*j)=xdist(ip,sign*(j-1));
00334           }
00335 
00336         //#if defined(DO_PROFILE)
00337         /*
00338           if (!underflow_flag)
00339           {
00340             hess_routine_and_constraint(ip,g,fg);  // calculate the hessian
00341                                             // at the conditional max
00342           }
00343 
00344           if (!underflow_flag)
00345           {
00346             ldet(ip,sign*j)=projected_hess_determinant(g,underflow_flag,
00347               xscale,ln_det_proj_jac(ip,sign*j));
00348           }
00349           else
00350           {
00351             ldet(ip,sign*j)=ldet(ip,sign*(j-1));
00352             ln_det_proj_jac(ip,sign*j)=ln_det_proj_jac(ip,sign*(j-1));
00353           }
00354           */
00355         //#endif
00356           }   // end of the j loop
00357         }
00358       }
00359       //for (int ix=-num_pp;ix<=num_pp;ix++)
00360       //{
00361        // xdist(ip,ix)=norm(elem_div(gprof(ix),xscale));
00362      // }
00363       if ( (option_match(ad_comm::argc,ad_comm::argv,"-prsave"))>-1)
00364       {
00365         adstring profrep_name=(likeprof_params::likeprofptr[ip]->label());
00366         ofstream ofs3((char*) (profrep_name+adstring(".pvl")));
00367         for (int ix=-num_pp;ix<=num_pp;ix++)
00368         {
00369           ofs3 << "#Step " << ix << endl;
00370           ofs3 << "#num sigmas " << all_num_sigs(ix) << endl;
00371           ofs3 << xmax(ix) << endl;
00372         }
00373       }
00374     }
00375   #if defined(DO_PROFILE)
00376     for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00377     {
00378       ldet(ip)=ldet(ip)-ln_det_proj_jac(ip);
00379     }
00380     {
00381       ofstream ofs("det.tmp");
00382       for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00383       {
00384         ofs << "the log dets" << endl;
00385         ofs << "ldet" << endl;
00386         ofs << ldet(ip) << endl << endl;
00387         ofs << "lndet_proj_jac" << endl;
00388         ofs << ln_det_proj_jac(ip) << endl << endl;
00389         ofs << "ldet-lndet_proj_jac" << endl;
00390         ofs << ldet(ip)-ln_det_proj_jac(ip) +ln_det_proj_jac(ip,0)
00391             << endl << endl;
00392       }
00393     }
00394     //for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00395     //{
00396      // ldet(ip)=ldet(ip)-(ln_det_proj_jac(ip)-ln_det_proj_jac(ip));
00397     //}
00398   #endif
00399       sigma=likeprof_params::likeprofptr[0]->get_sigma(); // this is the
00400       if (sigma==0.0)
00401         cerr << "Error standard dev of likeprof parameter is 0" << endl;
00402 #ifndef CURVE_CORRECT
00403     normalize_posterior_distribution(udet,siglevel,ofs2,num_pp,
00404       all_values,actual_value,global_min,offset,lprof,ldet,xdist,
00405       penalties);
00406 #else
00407     normalize_posterior_distribution(udet,siglevel,ofs2,num_pp,
00408       all_values,actual_value,global_min,offset,lprof,ldet,xdist,
00409       eigenvals,curvcor);
00410 #endif
00411   }
00412 
00413 /*
00414 ofstream of5("eigv.rpt");
00415 void get_ee(const dmatrix& hh, const ofstream& _of5)
00416 {
00417   ofstream& of5= (ofstream&) _of5;
00418   int mmin=hh.rowmin();
00419   int mmax=hh.rowmax();
00420   dvector l(mmin,mmax);
00421   dmatrix tmp(mmin,mmax,1,2);
00422   dvector ll(mmin,mmax);
00423   dmatrix e=eigenvectors(hh,l);
00424   for (int i=mmin;i<=mmax;i++)
00425   {
00426     ll(i)=e(i)*(hh*e(i));
00427     of5 << l(i) << "  " << ll(i) << endl;
00428   }
00429 }
00430 */