ADMB Documentation  11.1.2503
 All Classes Files Functions Variables Typedefs Friends Defines
mod_pmin.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_pmin.cpp 2432 2014-09-30 21:39:05Z 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     int ii=1;
00112     initial_params::copy_all_values(xsave,ii);
00113     double old_value;  // this is where we were
00114     double new_value;  // this is where we want to go
00115     double fprof = 0.0;
00116     double current_value;
00117     dmatrix lprof(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00118     dmatrix ldet(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00119     dmatrix ln_det_proj_jac(0,likeprof_params::num_likeprof_params-1,
00120       -num_pp,num_pp);
00121     dmatrix pdet(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00122     dmatrix actual_value(0,likeprof_params::num_likeprof_params-1,
00123       -num_pp,num_pp);
00124     dvector all_values(-num_pp,num_pp);
00125     dvector all_num_sigs(-num_pp,num_pp);
00126     //dvector xxxtmp(-num_pp,num_pp);
00127     //d3_array hesses(-num_pp,num_pp,1,nvar,1,nvar);
00128     dmatrix lg_jacob(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00129     dmatrix lg_prjacob(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00130     //dmatrix xxtmp(-num_pp,num_pp,1,nvar);
00131     dvector xvector(1,nvar);
00132     dmatrix xmax(-num_pp,num_pp,1,nvar);  // this holds the conditional max
00133     dmatrix gprof(-num_pp,num_pp,1,nvar);  // this holds the conditional max
00134     dmatrix fgrads(-num_pp,num_pp,1,nvar);  // this holds the conditional max
00135     dmatrix xdist(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00136     int sign=0;
00137     double sigma;
00138     ofstream ofs2((char*) (ad_comm::adprogram_name + adstring(".prf")) );
00139     double udet =unrestricted_hess_determinant();
00140     const int offset=0;
00141     dvector xscale(1,nvar);   // need to get scale from somewhere
00142     dvector likeprof_save(0,likeprof_params::num_likeprof_params-1);
00143     dmatrix penalties(0,likeprof_params::num_likeprof_params-1,-num_pp,num_pp);
00144     penalties.initialize();
00145 #ifdef CURVE_CORRECT
00146     int nlp=likeprof_params::num_likeprof_params;
00147     d3_array eigenvals(0,nlp-1,-num_pp,num_pp,1,nvar-1);
00148     d3_array curvcor(0,nlp-1,-num_pp,num_pp,1,nvar-1);
00149 #endif
00150 
00151     int ip;
00152     for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00153     {
00154       likeprof_save(ip)=likeprof_params::likeprofptr[ip]->get_value();
00155     }
00156     double final_weight;
00157     for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00158     {
00159       int sno=likeprof_params::likeprofptr[ip]->get_stepnumber();
00160       double snz=likeprof_params::likeprofptr[ip]->get_stepsize();
00161       if (sno)
00162       {
00163         num_pp=sno;
00164       }
00165       //snz (step size) must be greater than zero.
00166       //Check "void likeprof_params::set_stepsize(double x)"
00167       const double relsig = snz > 0 ? snz : 0.5;
00168 
00169       if (ip>0)
00170       {
00171         int ii=1;
00172         initial_params::restore_all_values(xsave,ii);
00173       }
00174       sigma=likeprof_params::likeprofptr[ip]->get_sigma(); // this is the
00175       if (sigma==0.0)
00176         cerr << "error standard dev of likeporf parameter is 0" << endl;
00177                                         // estimated sd
00178       old_value=likeprof_save(ip);
00179       old_value=old_value+offset*relsig*sigma;  // this is where we
00180       int bigint_flag=0;
00181       int bigint_flag1=0;
00182       double ldiff=0.0;
00183       double num_sigs;
00184       for (int i=1;i<=2;i++)  // go in positive and negative directions
00185       {
00186         num_sigs=0.0;
00187         bigint_flag=0;
00188         bigint_flag1=0;
00189         int underflow_flag=0;
00190         if (i>1) // get the parameter values at the global minimum
00191         {
00192           int ii=1;
00193           initial_params::restore_all_values(xsave,ii);
00194         }
00195         if (i==1)
00196         {
00197           sign=1;
00198         }
00199         else
00200         {
00201           sign=-1;
00202         }
00203         current_value=old_value;  // initialize at the minimum
00204         for (int j=0;j<=num_pp;j++)  // go in positive and negative directions
00205         {
00206           if (j!=0 || sign > 0)
00207           {
00208             if (bigint_flag==0)
00209             {
00210               num_sigs+=mult_factor(j)*relsig*sign;
00211               current_value+=mult_factor(j)*relsig*sign*sigma;
00212               new_value=current_value;
00213               // new_value=current_value+j*relsig*sign*sigma;
00214             }
00215             else
00216             {
00217               if (bigint_flag1==0)
00218               {
00219                 num_sigs+=1.5*relsig*sign;
00220                 current_value+=1.5*relsig*sign*sigma;
00221                 new_value=current_value;
00222               }
00223               else
00224               {
00225                 num_sigs+=2.5*relsig*sign;
00226                 current_value+=2.5*relsig*sign*sigma;
00227                 new_value=current_value;
00228               }
00229             }
00230 #if defined(USE_ADPVM)
00231           if (ad_comm::pvm_manager)
00232           {
00233             if (ad_comm::pvm_manager->mode==1) // master
00234             {
00235               send_int_to_slaves(3);
00236               pvm_master_prof_minimize(ip,sigma,new_value,fprof,underflow_flag,
00237                 global_min,penalties(ip,j*sign),final_weight); // get the
00238             }
00239           }
00240           else
00241 #endif
00242           {
00243             if (random_effects_flag==0)
00244             {
00245               prof_minimize(ip,sigma,new_value,fprof,underflow_flag,
00246               global_min,penalties(ip,j*sign),final_weight); // get the
00247                                                         // conditional max
00248             }
00249             else
00250             {
00251               prof_minimize_re(ip,sigma,new_value,fprof,underflow_flag,
00252               global_min,penalties(ip,j*sign),final_weight); // get the
00253                                                         // conditional max
00254             }
00255           }
00256           all_num_sigs(j*sign)=num_sigs;
00257           initial_params::xinit(xvector); // save the
00258           int ic=1;
00259           // save the conditional maximum
00260           initial_params::copy_all_values(xmax(sign*j),ic);
00261           /*int check=*/initial_params::stddev_scale(xscale,xvector);
00262         //#if defined(DO_PROFILE)
00263           //dvector curvscale(1,nvar);   // need to get scale from somewhere
00264         //#endif
00265 
00266        // #if defined(DO_PROFILE)
00267           // check=initial_params::stddev_curvscale(curvscale,xvector);
00268        // #endif
00269           //cout << "xscale = " << endl << xscale << endl;
00270           //{
00271            // ofstream ofs("xscale");
00272            // ofs << xscale << endl;
00273          // }
00274 #if defined(USE_DDOUBLE)
00275           lg_jacob(ip,sign*j)=sum(log(fabs(xscale)+double(1.e-60)));
00276 #else
00277           lg_jacob(ip,sign*j)=sum(log(fabs(xscale)+1.e-60));
00278 #endif
00279           if (!underflow_flag)
00280           {
00281             lprof(ip,sign*j)=fprof;
00282           }
00283           else
00284           {
00285             lprof(ip,sign*j)=lprof(ip,sign*(j-1))+2;
00286           }
00287           double xx=likeprof_params::likeprofptr[ip]->get_value();
00288           if (!underflow_flag)
00289           {
00290             actual_value(ip,sign*j)=xx;
00291           }
00292           else
00293           {
00294             actual_value(ip,sign*j)=new_value;  // this is where we
00295           }
00296 
00297           ldiff=fprof-lprof(ip,0);
00298           if ( ldiff > 40.0)
00299           {
00300             underflow_flag=1;
00301           }
00302 
00303           //if ( (fprof-lprof(ip,0)) > 3.0 && bigint_flag ==0)
00304           if (abs(j) > 4)
00305           {
00306             bigint_flag=1;
00307           }
00308           if (abs(j) > 5)
00309           {
00310             bigint_flag1=1;
00311           }
00312           // get the gradient for the profile likelihood variable at
00313           // the conditional maximum
00314           initial_params::current_phase = initial_params::max_number_phases;
00315           int nvar=initial_params::nvarcalc();
00316           dvector g(1,nvar);
00317           dvector fg(1,nvar);
00318           if (!underflow_flag)
00319           {
00320             // g is the grad. wrt the prof lik var
00321             // fg is the grad. wrt the log-likelihood
00322             get_particular_grad(ip,nvar,fg,g);
00323             gprof(sign*j)=g;
00324 
00325             fgrads(sign*j)=fg;
00326             xdist(ip,sign*j)=norm(elem_div(gprof(sign*j),xscale));
00327           }
00328           else
00329           {
00330             gprof(sign*j)=gprof(sign*(j-1));
00331             xdist(ip,sign*j)=xdist(ip,sign*(j-1));
00332           }
00333 
00334         //#if defined(DO_PROFILE)
00335         /*
00336           if (!underflow_flag)
00337           {
00338             hess_routine_and_constraint(ip,g,fg);  // calculate the hessian
00339                                             // at the conditional max
00340           }
00341 
00342           if (!underflow_flag)
00343           {
00344             ldet(ip,sign*j)=projected_hess_determinant(g,underflow_flag,
00345               xscale,ln_det_proj_jac(ip,sign*j));
00346           }
00347           else
00348           {
00349             ldet(ip,sign*j)=ldet(ip,sign*(j-1));
00350             ln_det_proj_jac(ip,sign*j)=ln_det_proj_jac(ip,sign*(j-1));
00351           }
00352           */
00353         //#endif
00354           }   // end of the j loop
00355         }
00356       }
00357       //for (int ix=-num_pp;ix<=num_pp;ix++)
00358       //{
00359        // xdist(ip,ix)=norm(elem_div(gprof(ix),xscale));
00360      // }
00361       if ( (option_match(ad_comm::argc,ad_comm::argv,"-prsave"))>-1)
00362       {
00363         adstring profrep_name=(likeprof_params::likeprofptr[ip]->label());
00364         ofstream ofs3((char*) (profrep_name+adstring(".pvl")));
00365         for (int ix=-num_pp;ix<=num_pp;ix++)
00366         {
00367           ofs3 << "#Step " << ix << endl;
00368           ofs3 << "#num sigmas " << all_num_sigs(ix) << endl;
00369           ofs3 << xmax(ix) << endl;
00370         }
00371       }
00372     }
00373   #if defined(DO_PROFILE)
00374     for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00375     {
00376       ldet(ip)=ldet(ip)-ln_det_proj_jac(ip);
00377     }
00378     {
00379       ofstream ofs("det.tmp");
00380       for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00381       {
00382         ofs << "the log dets" << endl;
00383         ofs << "ldet" << endl;
00384         ofs << ldet(ip) << endl << endl;
00385         ofs << "lndet_proj_jac" << endl;
00386         ofs << ln_det_proj_jac(ip) << endl << endl;
00387         ofs << "ldet-lndet_proj_jac" << endl;
00388         ofs << ldet(ip)-ln_det_proj_jac(ip) +ln_det_proj_jac(ip,0)
00389             << endl << endl;
00390       }
00391     }
00392     //for (ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00393     //{
00394      // ldet(ip)=ldet(ip)-(ln_det_proj_jac(ip)-ln_det_proj_jac(ip));
00395     //}
00396   #endif
00397       sigma=likeprof_params::likeprofptr[0]->get_sigma(); // this is the
00398       if (sigma==0.0)
00399         cerr << "Error standard dev of likeprof parameter is 0" << endl;
00400 #ifndef CURVE_CORRECT
00401     normalize_posterior_distribution(udet,siglevel,ofs2,num_pp,
00402       all_values,actual_value,global_min,offset,lprof,ldet,xdist,
00403       penalties);
00404 #else
00405     normalize_posterior_distribution(udet,siglevel,ofs2,num_pp,
00406       all_values,actual_value,global_min,offset,lprof,ldet,xdist,
00407       eigenvals,curvcor);
00408 #endif
00409   }
00410 
00411 /*
00412 ofstream of5("eigv.rpt");
00413 void get_ee(const dmatrix& hh, const ofstream& _of5)
00414 {
00415   ofstream& of5= (ofstream&) _of5;
00416   int mmin=hh.rowmin();
00417   int mmax=hh.rowmax();
00418   dvector l(mmin,mmax);
00419   dmatrix tmp(mmin,mmax,1,2);
00420   dvector ll(mmin,mmax);
00421   dmatrix e=eigenvectors(hh,l);
00422   for (int i=mmin;i<=mmax;i++)
00423   {
00424     ll(i)=e(i)*(hh*e(i));
00425     of5 << l(i) << "  " << ll(i) << endl;
00426   }
00427 }
00428 */