ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
newmodmn.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: newmodmn.cpp 2587 2014-11-07 19:38:15Z 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 polint(double * xa,double * ya,double x);
00010 
00011 void get_confidence_interval(const dvector& left_bd, const dvector& right_bd,
00012   dmatrix& ms, const dvector& xs, const dvector& siglevel,
00013   const int& level_index, int index);
00014 void get_onesided_intervals(const dvector& left_bd, const dvector& right_bd,
00015   dmatrix& ms, const dvector& xs, const dvector& siglevel,
00016   const int& level_index, int index);
00017 void report_confidence_limits(const ofstream& ofs3,int numsig_levels,
00018   const dvector& siglevel, const dvector& left_bd, const dvector& right_bd);
00019 void report_onesided_confidence_limits(const ofstream& ofs3,int numsig_levels,
00020   const dvector& siglevel, const dvector& left_bd, const dvector& right_bd,
00021   int ip);
00022 
00023 void report_confidence_limits(const ofstream& _ofs3,int numsig_levels,
00024   const dvector& siglevel, const dvector& left_bd, const dvector& right_bd)
00025 {
00026   ofstream& ofs3=(ofstream&) _ofs3;
00027   ofs3 << "Minimum width confidence limits:" << endl
00028        <<                    " significance level  lower bound"
00029        << "  upper bound" << endl;
00030   for (int i=1;i<=numsig_levels;i++)
00031   {
00032     ofs3 << "                 " << setw(12) << siglevel(i)
00033          << "             " << left_bd(i) << "     " << right_bd(i) << endl;
00034   }
00035   ofs3 << endl;
00036 }
00037 
00038 void report_onesided_confidence_limits(const ofstream& _ofs3,int numsig_levels,
00039   const dvector& siglevel, const dvector& left_bd, const dvector& right_bd,
00040   int ip)
00041 {
00042   ofstream& ofs3=(ofstream&) _ofs3;
00043   int i;
00044   for (i=1;i<=numsig_levels;i++)
00045   {
00046     ofs3 << "The probability is " << setw(7) << siglevel(i) << " that "
00047          << likeprof_params::likeprofptr[ip]->label()
00048          << " is greater than " << left_bd(i) << endl;
00049   }
00050   ofs3 << endl;
00051   for (i=1;i<=numsig_levels;i++)
00052   {
00053     ofs3 << "The probability is " << setw(7) << siglevel(i) << " that "
00054          << likeprof_params::likeprofptr[ip]->label()
00055          << " is less than " << right_bd(i) << endl;
00056   }
00057   ofs3 << endl;
00058 }
00059 
00060   dvector smooth(const dvector& v);
00061 
00062 #ifndef CURVE_CORRECT
00063   void function_minimizer::normalize_posterior_distribution(double udet,
00064     const dvector& siglevel, const ofstream& _ofs2,int num_pp,
00065     const dvector& _all_values, const dmatrix& actual_value,double global_min,
00066     int offset, const dmatrix& lprof, const dmatrix& ldet, const dmatrix& xdist,
00067     const dmatrix& penalties)
00068   /*
00069   void function_minimizer::normalize_posterior_distribution(double udet,
00070     dvector& siglevel, const ofstream& ofs2,int num_pp, const dvector& all_values,
00071     dmatrix& actual_value,double global_min,int offset, const dmatrix& lprof,
00072     dmatrix& ldet, const dmatrix& xdist, const dmatrix& penalties)
00073    */
00074    // dmatrix& ldet, const dmatrix& xdist, const dmatrix& penalties,
00075    // const dmatrix& lg_jacob)
00076 #else
00077 
00078   void function_minimizer::normalize_posterior_distribution(double udet,
00079     const dvector& siglevel, const ofstream& ofs2,int num_pp,
00080     const dvector& all_values, const dmatrix& actual_value,
00081     double global_min,
00082     int offset, const dmatrix& lprof, const dmatrix& ldet,
00083     const dmatrix& xdist,
00084     const d3_array& eigenvals,d3_array& curvcor)
00085 /*
00086   void function_minimizer::normalize_posterior_distribution(double udet,
00087     const dvector& siglevel, const ofstream& ofs2,int num_pp,
00088     const dvector& all_values,
00089     const dmatrix& actual_value,double global_min,int offset,
00090     const dmatrix& lprof,
00091     const dmatrix& ldet, const dmatrix& xdist,
00092     const d3_array& eigenvals,d3_array& curvcor)
00093 */
00094 #endif
00095   {
00096     dvector& all_values=(dvector&) _all_values;
00097     ofstream& ofs2=(ofstream&) _ofs2;
00098     ofstream savef("diags");
00099     //ofstream ofs3((char*) (ad_comm::adprogram_name + ".plt") );
00100     //dvector siglevel("{.90,.95,.975}");
00101 #ifndef CURVE_CORRECT
00102     dmatrix left_bd(1,3,1,3);
00103     dmatrix right_bd(1,3,1,3);
00104     dmatrix lower_bd(1,3,1,3);
00105     dmatrix upper_bd(1,3,1,3);
00106 #else
00107     dmatrix left_bd(0,3,1,3);
00108     dmatrix right_bd(0,3,1,3);
00109     dmatrix lower_bd(0,3,1,3);
00110     dmatrix upper_bd(0,3,1,3);
00111 #endif
00112     double sigma;
00113     /*int nvar=*/initial_params::nvarcalc();
00114     int numsig_levels=siglevel.indexmax();
00115     for (int ip=0;ip<likeprof_params::num_likeprof_params;ip++)
00116     {
00117      {
00118       adstring profrep_name=(likeprof_params::likeprofptr[ip]->label());
00119       size_t llen = length(profrep_name);
00120       if (llen>8) profrep_name=profrep_name(1,8);
00121       ofstream ofs3((char*) (profrep_name+adstring(".plt")));
00122       sigma=likeprof_params::likeprofptr[ip]->get_sigma();
00123       //double diff;
00124       ofs2 << likeprof_params::likeprofptr[ip]->label() << ":" << endl;
00125       ofs3 << likeprof_params::likeprofptr[ip]->label() << ":" << endl;
00126       ofs2 << sigma << endl  << global_min << " " << udet << endl;
00127       dvector tempint0(-num_pp,num_pp);
00128       dvector tempint1(-num_pp,num_pp);
00129       dvector tempint2(-num_pp,num_pp);
00130       {
00131         ofstream ofs("dgs2");
00132         ofs << "lprof" << endl << lprof << endl;
00133       }
00134       for (int j=-num_pp;j<=num_pp;j++) //go in positive and negative directions
00135       {
00136         all_values(j)=actual_value(ip,j-offset);
00137         double lp=lprof(ip,j);
00138       #if defined(DO_PROFILE)
00139         tempint0(j)= exp(global_min-lp+.5*(-ldet(ip,j)+ldet(ip,0)));
00140       #endif
00141         tempint1(j)= exp(global_min-lp);
00142 
00143         tempint2(j)=
00144           square((actual_value(ip,j)-actual_value(ip,-offset))/
00145             (1.414*sigma));
00146         tempint2(j)=exp(-tempint2(j));
00147       }
00148       dmatrix m(1,3,-num_pp,num_pp);
00149       for (int j=-num_pp;j<=num_pp;j++)
00150       {
00151       #if defined(DO_PROFILE)
00152         m(1,j)=tempint0(j)/xdist(ip,j);
00153       #endif
00154         if (xdist(ip,j)<1.e-50)
00155         {
00156           cerr << " xdist(" << ip << "," << j << ") too small" << endl;
00157           char ch;
00158           cin >> ch;
00159           m(2,j)=1.e+20;
00160         }
00161         else
00162         {
00163           // profile likelihood adjusted for gradient of dep var
00164           m(2,j)=tempint1(j)/xdist(ip,j);
00165         }
00166         //m(2,j)=tempint1(j);
00167       }
00168       //m(2,num_pp)=tempint1(num_pp)/xdist(ip,num_pp);
00169       m(3)=tempint2; //+ 1.e-4*max(tempint2);
00170 
00171      /*
00172       savef << "penalties" << endl << setw(9) << setprecision(4)
00173             << penalties(ip) << endl;
00174       savef << "normalized exp(lg_jacob)" << endl << setw(9) << setprecision(4)
00175             << exp(2.*(lg_jacob(ip)-lg_jacob(ip,0))) << endl;
00176      */
00177       savef << "tempint1 " << endl << setw(9) << setprecision(3)
00178             << tempint1 << endl;
00179      #if defined(DO_PROFILE)
00180       savef << "m(1) " << endl << setw(9) << setprecision(3) << m(1) << endl;
00181      #endif
00182       savef << "m(2) " << endl << setw(9) << setprecision(3) << m(2) << endl;
00183       savef << "m(3) " << endl << setw(9) << setprecision(3) << m(3) << endl;
00184       savef << "xdistance" << endl;
00185       savef << xdist(ip) << endl << endl;
00186 
00187       int min_ind = -num_pp;
00188       int max_ind = num_pp;
00189 
00190       dvector xs(7*min_ind,7*max_ind);
00191       dmatrix ms(1,3,7*min_ind,7*max_ind);
00192 
00193       int kmult=7;
00194       for (int k=min_ind;k<=max_ind;k++)
00195       {
00196         double val=all_values(k);
00197         xs(kmult*k)=val;
00198         if (k<max_ind)
00199         {
00200           double diff=all_values(k+1)-all_values(k);
00201           for (int i=1;i<kmult;i++)
00202           {
00203             xs(kmult*k+i)=val+i/double(kmult)*diff;
00204           }
00205         }
00206       }
00207       {
00208         int mmin=m.rowmin();
00209         int mmax=m.rowmax();
00210         for (int i=mmin;i<=mmax;i++)
00211         {
00212           int cmin=m(i).indexmin();
00213           int cmax=m(i).indexmax();
00214           for (int j=cmin;j<=cmax;j++)
00215           {
00216             if (m(i,j)<=0.0) m(i,j)=1.e-50;
00217           }
00218         }
00219       }
00220       //dmatrix lm=log(m);
00221       dmatrix lm=m;
00222       int lowlimit=2;
00223       #if defined(DO_PROFILE)
00224         lowlimit=1;
00225       #else
00226         lowlimit=2;
00227       #endif
00228 
00229       for (int ii=lowlimit;ii<=3;ii++)
00230       {
00231         int k;
00232         for (k=min_ind;k<=0;k++)
00233         {
00234           //double * xa=(&all_values(k))-1;
00235           //double * ya=(&lm(ii,k))-1;
00236           //ms(ii,kmult*k)=exp(lm(ii,k));
00237           ms(ii,kmult*k)=lm(ii,k);
00238           if (k<max_ind)
00239           {
00240             double l=lm(ii,k);
00241             double u=lm(ii,k+1);
00242             for (int i=1;i<kmult;i++)
00243             {
00244               ms(ii,kmult*k+i)=l+double(i)/kmult*(u-l);
00245               //ms(ii,kmult*k+i)=polint(xa,ya,xs(k*kmult+i));
00246             }
00247           }
00248         }
00249 
00250         for (k=1;k<=max_ind;k++)
00251         {
00252           //double * xa=(&all_values(k))-2;
00253           //double * ya=(&lm(ii,k))-2;
00254           //ms(ii,kmult*k)=exp(lm(ii,k));
00255           ms(ii,kmult*k)=lm(ii,k);
00256           if (k<max_ind)
00257           {
00258             double l=lm(ii,k);
00259             double u=lm(ii,k+1);
00260             for (int i=1;i<kmult;i++)
00261             {
00262               ms(ii,kmult*k+i)=l+double(i)/kmult*(u-l);
00263             }
00264           }
00265         }
00266       }
00267 
00268     /*
00269       savef << "ms(2) " << endl << setw(9) << setprecision(3) << ms(2) << endl;
00270       savef << "ms(3) " << endl << setw(9) << setprecision(3) << ms(3) << endl;
00271     */
00272 
00273       dvector ssum(1,3);
00274       ssum.initialize();
00275       for (int j=lowlimit;j<=3;j++)
00276       {
00277         for (int i=7*min_ind;i<7*max_ind;i++)
00278         {
00279           if (ms(j,i)<0.0)
00280           {
00281             ms(j,i)=0.0;
00282           }
00283           else
00284           {
00285             ssum(j)+=ms(j,i)*(xs(i+1)-xs(i));
00286           }
00287         }
00288       }
00289       for (int j=lowlimit;j<=3;j++)
00290       {
00291         if (ssum(j) !=0)
00292         {
00293         /*
00294           cout << ms(j) << endl << ssum(j) << endl << endl;
00295           char ch;
00296           cin >> ch;
00297         */
00298           ms(j)/=ssum(j);
00299         }
00300       }
00301       // now do the integrals
00302       for (int j=lowlimit;j<=3;j++)
00303       {
00304         int level_index=1;
00305         do
00306         {
00307           get_confidence_interval(left_bd(j),right_bd(j),ms,xs,
00308             siglevel,level_index,j);
00309           get_onesided_intervals(lower_bd(j),upper_bd(j),ms,xs,
00310             siglevel,level_index,j);
00311           level_index++;
00312         }
00313         while (level_index <= numsig_levels);
00314       }
00315 
00316       double min_cutoff = 1.e-3/sigma;
00317       int i;
00318       for (i=7*min_ind;i<=0;i++)
00319       {
00320         if (max(ms(2,i),ms(3,i)) > min_cutoff) break;
00321       }
00322       int new_min_ind = int(max(i,7*min_ind));
00323       for (i=0;i<=7*max_ind;i++)
00324       {
00325         if (max(ms(2,i),ms(3,i)) < min_cutoff) break;
00326       }
00327       int new_max_ind = min(i,7*max_ind);
00328       dmatrix output(1,2,new_min_ind,new_max_ind);
00329 
00330       output(1)=xs(new_min_ind,new_max_ind);
00331 
00332     #if defined(DO_PROFILE)
00333       output(2)=ms(1)(new_min_ind,new_max_ind);
00334       {
00335         ofs3 << "Adjusted Profile likelihood" << endl;
00336         ofs3 << trans(output) << endl;
00337         report_confidence_limits(ofs3,numsig_levels,siglevel,left_bd(1),
00338           right_bd(1));
00339         ofs3 << "One sided confidence limits for the "
00340                 "adjusted profile likelihood:" << endl << endl;
00341         report_onesided_confidence_limits(ofs3,numsig_levels,siglevel,
00342           lower_bd(1),upper_bd(1),ip);
00343       }
00344      #endif
00345 
00346       output(2)=ms(2)(new_min_ind,new_max_ind);
00347       {
00348         ofs3 << "Profile likelihood" << endl;
00349         ofs3 << trans(output) << endl;
00350         report_confidence_limits(ofs3,numsig_levels,siglevel,left_bd(2),
00351           right_bd(2));
00352         ofs3 << "One sided confidence limits for the "
00353                 "profile likelihood:" << endl << endl;
00354         report_onesided_confidence_limits(ofs3,numsig_levels,siglevel,
00355           lower_bd(2),upper_bd(2),ip);
00356       }
00357 
00358       output(2)=ms(3)(new_min_ind,new_max_ind);
00359       {
00360         ofs3 << "Normal approximation" << endl;
00361         ofs3 << trans(output) << endl;
00362         report_confidence_limits(ofs3,numsig_levels,siglevel,left_bd(3),
00363           right_bd(3));
00364         ofs3 << "One sided confidence limits for the "
00365                 "Normal approximation:" << endl << endl;
00366         report_onesided_confidence_limits(ofs3,numsig_levels,siglevel,
00367           lower_bd(3),upper_bd(3),ip);
00368       }
00369      }
00370     }
00371   }
00372 
00373 
00374   dvector smooth(const dvector& v)
00375   {
00376     int mmin=v.indexmin();
00377     int mmax=v.indexmax();
00378     int diff=mmax-mmin+1;
00379     dvector tmp(mmin,mmax);
00380     tmp(mmin)=.8*v(mmin)+.2*v(mmin+1);
00381     tmp(mmin+1)=.2*v(mmin)+.6*v(mmin+1)+.2*v(mmin+2);
00382     for (int i=mmin+2;i<=mmin+diff/4;i++)
00383     {
00384       tmp(i)=.1*v(i-2)+.2*v(i-1)+.4*v(i)+.2*v(i+1)+.1*v(i+2);
00385     }
00386     for (int i=mmin+diff/4+1;i<mmax-diff/4;i++)
00387     {
00388       tmp(i)=v(i);
00389     }
00390     for (int i=mmax-diff/4;i<=mmax-2;i++)
00391     {
00392       tmp(i)=.1*v(i-2)+.2*v(i-1)+.4*v(i)+.2*v(i+1)+.1*v(i+2);
00393     }
00394     tmp(mmax)=.8*v(mmax)+.2*v(mmax-1);
00395     tmp(mmax-1)=.2*v(mmax)+.6*v(mmax-1)+.2*v(mmax-2);
00396     return tmp;
00397   }
00398 
00399 double polint(double * xa,double * ya,double x)
00400 {
00401   double prod1=(xa[1]-xa[2])*(xa[1]-xa[3]);
00402   double prod2=(xa[2]-xa[1])*(xa[2]-xa[3]);
00403   double prod3=(xa[3]-xa[1])*(xa[3]-xa[2]);
00404   if (prod1==0)
00405   {
00406     double y=ya[1];
00407     return y;
00408   }
00409   if (prod2==0)
00410   {
00411     double y=ya[2];
00412     return y;
00413   }
00414   if (prod3==0)
00415   {
00416     double y=ya[2];
00417     return y;
00418   }
00419   double y= (x-xa[2])*(x-xa[3])/prod1*ya[1]
00420            +(x-xa[1])*(x-xa[3])/prod2*ya[2]
00421            +(x-xa[1])*(x-xa[2])/prod3*ya[3];
00422   return y;
00423 }