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