ADMB Documentation  11.1.2497
 All Classes Files Functions Variables Typedefs Friends Defines
newmodmn.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: newmodmn.cpp 1919 2014-04-22 22:02:01Z 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       int 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       int j;
00135       for (j=-num_pp;j<=num_pp;j++) //go in positive and negative directions
00136       {
00137         all_values(j)=actual_value(ip,j-offset);
00138         double lp=lprof(ip,j);
00139       #if defined(DO_PROFILE)
00140         tempint0(j)= exp(global_min-lp+.5*(-ldet(ip,j)+ldet(ip,0)));
00141       #endif
00142         tempint1(j)= exp(global_min-lp);
00143 
00144         tempint2(j)=
00145           square((actual_value(ip,j)-actual_value(ip,-offset))/
00146             (1.414*sigma));
00147         tempint2(j)=exp(-tempint2(j));
00148       }
00149       dmatrix m(1,3,-num_pp,num_pp);
00150       for (j=-num_pp;j<=num_pp;j++)
00151       {
00152       #if defined(DO_PROFILE)
00153         m(1,j)=tempint0(j)/xdist(ip,j);
00154       #endif
00155         if (xdist(ip,j)<1.e-50)
00156         {
00157           cerr << " xdist(" << ip << "," << j << ") too small" << endl;
00158           char ch;
00159           cin >> ch;
00160           m(2,j)=1.e+20;
00161         }
00162         else
00163         {
00164           // profile likelihood adjusted for gradient of dep var
00165           m(2,j)=tempint1(j)/xdist(ip,j);
00166         }
00167         //m(2,j)=tempint1(j);
00168       }
00169       //m(2,num_pp)=tempint1(num_pp)/xdist(ip,num_pp);
00170       m(3)=tempint2; //+ 1.e-4*max(tempint2);
00171 
00172      /*
00173       savef << "penalties" << endl << setw(9) << setprecision(4)
00174             << penalties(ip) << endl;
00175       savef << "normalized exp(lg_jacob)" << endl << setw(9) << setprecision(4)
00176             << exp(2.*(lg_jacob(ip)-lg_jacob(ip,0))) << endl;
00177      */
00178       savef << "tempint1 " << endl << setw(9) << setprecision(3)
00179             << tempint1 << endl;
00180      #if defined(DO_PROFILE)
00181       savef << "m(1) " << endl << setw(9) << setprecision(3) << m(1) << endl;
00182      #endif
00183       savef << "m(2) " << endl << setw(9) << setprecision(3) << m(2) << endl;
00184       savef << "m(3) " << endl << setw(9) << setprecision(3) << m(3) << endl;
00185       savef << "xdistance" << endl;
00186       savef << xdist(ip) << endl << endl;
00187 
00188       int min_ind = -num_pp;
00189       int max_ind = num_pp;
00190 
00191       dvector xs(7*min_ind,7*max_ind);
00192       dmatrix ms(1,3,7*min_ind,7*max_ind);
00193 
00194       int kmult=7;
00195       for (int k=min_ind;k<=max_ind;k++)
00196       {
00197         double val=all_values(k);
00198         xs(kmult*k)=val;
00199         if (k<max_ind)
00200         {
00201           double diff=all_values(k+1)-all_values(k);
00202           for (int i=1;i<kmult;i++)
00203           {
00204             xs(kmult*k+i)=val+i/double(kmult)*diff;
00205           }
00206         }
00207       }
00208       {
00209         int mmin=m.rowmin();
00210         int mmax=m.rowmax();
00211         for (int i=mmin;i<=mmax;i++)
00212         {
00213           int cmin=m(i).indexmin();
00214           int cmax=m(i).indexmax();
00215           for (int j=cmin;j<=cmax;j++)
00216           {
00217             if (m(i,j)<=0.0) m(i,j)=1.e-50;
00218           }
00219         }
00220       }
00221       //dmatrix lm=log(m);
00222       dmatrix lm=m;
00223       int lowlimit=2;
00224       #if defined(DO_PROFILE)
00225         lowlimit=1;
00226       #else
00227         lowlimit=2;
00228       #endif
00229 
00230       for (int ii=lowlimit;ii<=3;ii++)
00231       {
00232         int k;
00233         for (k=min_ind;k<=0;k++)
00234         {
00235           //double * xa=(&all_values(k))-1;
00236           //double * ya=(&lm(ii,k))-1;
00237           //ms(ii,kmult*k)=exp(lm(ii,k));
00238           ms(ii,kmult*k)=lm(ii,k);
00239           if (k<max_ind)
00240           {
00241             double l=lm(ii,k);
00242             double u=lm(ii,k+1);
00243             for (int i=1;i<kmult;i++)
00244             {
00245               ms(ii,kmult*k+i)=l+double(i)/kmult*(u-l);
00246               //ms(ii,kmult*k+i)=polint(xa,ya,xs(k*kmult+i));
00247             }
00248           }
00249         }
00250 
00251         for (k=1;k<=max_ind;k++)
00252         {
00253           //double * xa=(&all_values(k))-2;
00254           //double * ya=(&lm(ii,k))-2;
00255           //ms(ii,kmult*k)=exp(lm(ii,k));
00256           ms(ii,kmult*k)=lm(ii,k);
00257           if (k<max_ind)
00258           {
00259             double l=lm(ii,k);
00260             double u=lm(ii,k+1);
00261             for (int i=1;i<kmult;i++)
00262             {
00263               ms(ii,kmult*k+i)=l+double(i)/kmult*(u-l);
00264             }
00265           }
00266         }
00267       }
00268 
00269     /*
00270       savef << "ms(2) " << endl << setw(9) << setprecision(3) << ms(2) << endl;
00271       savef << "ms(3) " << endl << setw(9) << setprecision(3) << ms(3) << endl;
00272     */
00273 
00274       dvector ssum(1,3);
00275       ssum.initialize();
00276       for (j=lowlimit;j<=3;j++)
00277       {
00278         for (int i=7*min_ind;i<7*max_ind;i++)
00279         {
00280           if (ms(j,i)<0.0)
00281           {
00282             ms(j,i)=0.0;
00283           }
00284           else
00285           {
00286             ssum(j)+=ms(j,i)*(xs(i+1)-xs(i));
00287           }
00288         }
00289       }
00290       for (j=lowlimit;j<=3;j++)
00291       {
00292         if (ssum(j) !=0)
00293         {
00294         /*
00295           cout << ms(j) << endl << ssum(j) << endl << endl;
00296           char ch;
00297           cin >> ch;
00298         */
00299           ms(j)/=ssum(j);
00300         }
00301       }
00302       // now do the integrals
00303       for (j=lowlimit;j<=3;j++)
00304       {
00305         int level_index=1;
00306         do
00307         {
00308           get_confidence_interval(left_bd(j),right_bd(j),ms,xs,
00309             siglevel,level_index,j);
00310           get_onesided_intervals(lower_bd(j),upper_bd(j),ms,xs,
00311             siglevel,level_index,j);
00312           level_index++;
00313         }
00314         while (level_index <= numsig_levels);
00315       }
00316 
00317       double min_cutoff = 1.e-3/sigma;
00318       int i;
00319       for (i=7*min_ind;i<=0;i++)
00320       {
00321         if (max(ms(2,i),ms(3,i)) > min_cutoff) break;
00322       }
00323       int new_min_ind = int(max(i,7*min_ind));
00324       for (i=0;i<=7*max_ind;i++)
00325       {
00326         if (max(ms(2,i),ms(3,i)) < min_cutoff) break;
00327       }
00328       int new_max_ind = min(i,7*max_ind);
00329       dmatrix output(1,2,new_min_ind,new_max_ind);
00330 
00331       output(1)=xs(new_min_ind,new_max_ind);
00332 
00333     #if defined(DO_PROFILE)
00334       output(2)=ms(1)(new_min_ind,new_max_ind);
00335       {
00336         ofs3 << "Adjusted Profile likelihood" << endl;
00337         ofs3 << trans(output) << endl;
00338         report_confidence_limits(ofs3,numsig_levels,siglevel,left_bd(1),
00339           right_bd(1));
00340         ofs3 << "One sided confidence limits for the "
00341                 "adjusted profile likelihood:" << endl << endl;
00342         report_onesided_confidence_limits(ofs3,numsig_levels,siglevel,
00343           lower_bd(1),upper_bd(1),ip);
00344       }
00345      #endif
00346 
00347       output(2)=ms(2)(new_min_ind,new_max_ind);
00348       {
00349         ofs3 << "Profile likelihood" << endl;
00350         ofs3 << trans(output) << endl;
00351         report_confidence_limits(ofs3,numsig_levels,siglevel,left_bd(2),
00352           right_bd(2));
00353         ofs3 << "One sided confidence limits for the "
00354                 "profile likelihood:" << endl << endl;
00355         report_onesided_confidence_limits(ofs3,numsig_levels,siglevel,
00356           lower_bd(2),upper_bd(2),ip);
00357       }
00358 
00359       output(2)=ms(3)(new_min_ind,new_max_ind);
00360       {
00361         ofs3 << "Normal approximation" << endl;
00362         ofs3 << trans(output) << endl;
00363         report_confidence_limits(ofs3,numsig_levels,siglevel,left_bd(3),
00364           right_bd(3));
00365         ofs3 << "One sided confidence limits for the "
00366                 "Normal approximation:" << endl << endl;
00367         report_onesided_confidence_limits(ofs3,numsig_levels,siglevel,
00368           lower_bd(3),upper_bd(3),ip);
00369       }
00370      }
00371     }
00372   }
00373 
00374 
00375   dvector smooth(const dvector& v)
00376   {
00377     int mmin=v.indexmin();
00378     int mmax=v.indexmax();
00379     int diff=mmax-mmin+1;
00380     dvector tmp(mmin,mmax);
00381     tmp(mmin)=.8*v(mmin)+.2*v(mmin+1);
00382     tmp(mmin+1)=.2*v(mmin)+.6*v(mmin+1)+.2*v(mmin+2);
00383     int i;
00384     for (i=mmin+2;i<=mmin+diff/4;i++)
00385     {
00386       tmp(i)=.1*v(i-2)+.2*v(i-1)+.4*v(i)+.2*v(i+1)+.1*v(i+2);
00387     }
00388     for (i=mmin+diff/4+1;i<mmax-diff/4;i++)
00389     {
00390       tmp(i)=v(i);
00391     }
00392     for (i=mmax-diff/4;i<=mmax-2;i++)
00393     {
00394       tmp(i)=.1*v(i-2)+.2*v(i-1)+.4*v(i)+.2*v(i+1)+.1*v(i+2);
00395     }
00396     tmp(mmax)=.8*v(mmax)+.2*v(mmax-1);
00397     tmp(mmax-1)=.2*v(mmax)+.6*v(mmax-1)+.2*v(mmax-2);
00398     return tmp;
00399   }
00400 
00401 double polint(double * xa,double * ya,double x)
00402 {
00403   double prod1=(xa[1]-xa[2])*(xa[1]-xa[3]);
00404   double prod2=(xa[2]-xa[1])*(xa[2]-xa[3]);
00405   double prod3=(xa[3]-xa[1])*(xa[3]-xa[2]);
00406   if (prod1==0)
00407   {
00408     double y=ya[1];
00409     return y;
00410   }
00411   if (prod2==0)
00412   {
00413     double y=ya[2];
00414     return y;
00415   }
00416   if (prod3==0)
00417   {
00418     double y=ya[2];
00419     return y;
00420   }
00421   double y= (x-xa[2])*(x-xa[3])/prod1*ya[1]
00422            +(x-xa[1])*(x-xa[3])/prod2*ya[2]
00423            +(x-xa[1])*(x-xa[2])/prod3*ya[3];
00424   return y;
00425 }