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