ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
newmodm3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: newmodm3.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 int xxxmax(int x,int y)
00010 {
00011   if (x<y) return y;
00012   return x;
00013 }
00014 int xxxmin(int x,int y)
00015 {
00016   if (x<y) return x;
00017   return y;
00018 }
00019 void get_confidence_interval(const dvector& _left_bd, const dvector& _right_bd,
00020   dmatrix& ms, const dvector& xs, const dvector& siglevel,
00021   const int& level_index, int index)
00022   {
00023     dvector& left_bd=(dvector&) _left_bd;
00024     dvector& right_bd=(dvector&) _right_bd;
00025     dvector& fval=ms(index);
00026     int lb=0;
00027     int rb=0;
00028     double xdiff=0.0;
00029     double isum=0;
00030 
00031     int mmin=fval.indexmin();
00032     int mmax=fval.indexmax();
00033     double tmp=fval(mmin);
00034     int imax=mmin;
00035     int i;
00036     for (i=mmin+1;i<=mmax;i++)
00037     {
00038       if (fval(i)>tmp)
00039       {
00040         tmp=fval(i);
00041         imax=i;
00042       }
00043     }
00044     if (imax>mmin)
00045     {
00046       lb=imax-1;
00047       rb=imax;
00048     }
00049     else
00050     {
00051       lb=imax;
00052       rb=imax+1;
00053     }
00054     double integral=0.0;
00055 
00056     for (i=mmin+1;i<=mmax;i++)
00057     {
00058       integral+=fval(i)*(xs(i)-xs(i-1));
00059     }
00060     cout << integral << endl;
00061     do
00062     {
00063       if (lb <= fval.indexmin() &&  rb > fval.indexmax())
00064       {
00065         int lower=xxxmax(lb,fval.indexmin());
00066         int upper=xxxmin(rb-1,fval.indexmax());
00067         left_bd(level_index)=xs(lower);
00068         right_bd(level_index)=xs(upper);
00069         break;
00070       }
00071 
00072       if (fval(lb)>=fval(rb) || rb==fval.indexmax())
00073       {
00074         if (lb>fval.indexmin())
00075         {
00076           xdiff=xs(lb)-xs(lb-1);
00077         }
00078         else
00079         {
00080           int minind=fval.indexmin();
00081           xdiff=xs(minind+1)-xs(minind);
00082         }
00083         double incr=fval(lb)*xdiff/integral;
00084         //double incr=fval(lb);
00085         if ( incr >= (siglevel(level_index)-isum))
00086         {
00087           double diff = siglevel(level_index) - isum;
00088           double delta=diff/incr;
00089           if (lb>fval.indexmin())
00090           {
00091             left_bd(level_index)=xs(lb)-delta*(xdiff);
00092           }
00093           else
00094           {
00095             left_bd(level_index)=xs(fval.indexmin());
00096           }
00097           right_bd(level_index)=xs(rb);
00098           break;
00099         }
00100         isum+=incr;
00101         lb--;
00102       }
00103       else
00104       {
00105         xdiff=xs(rb)-xs(rb-1);
00106         double incr=fval(rb)*xdiff/integral;
00107         //double incr=fval(rb);
00108         if ( incr >= (siglevel(level_index)-isum) ||
00109           rb == mmax)
00110         {
00111           double diff = siglevel(level_index) - isum;
00112           double delta=diff/incr;
00113           left_bd(level_index)=xs(lb);
00114           right_bd(level_index)=xs(rb)+delta*(xdiff);
00115           break;
00116         }
00117         isum+=incr;
00118         rb++;
00119       }
00120     }
00121     while (1);
00122   }
00123 
00124 void get_onesided_intervals(const dvector& _left_bd, const dvector& _right_bd,
00125   dmatrix& ms, const dvector& xs, const dvector& siglevel,
00126   const int& level_index, int index)
00127   {
00128     dvector& left_bd=(dvector&) _left_bd;
00129     dvector& right_bd=(dvector&) _right_bd;
00130     dvector& fval=ms(index);
00131     int lb=fval.indexmin()+1;
00132     int ub=fval.indexmax();
00133     double xdiff=0.0;
00134     double isum=0;
00135     double tmpsum=0.0;
00136     int ii;
00137     for (ii=lb+1;ii<=ub;ii++)
00138     {
00139       tmpsum+=fval(ii)*(xs(ii)-xs(ii-1));
00140     }
00141 
00142     isum=0.0;
00143     for (ii=lb+1;ii<=ub;ii++)
00144     {
00145       xdiff=xs(ii)-xs(ii-1);
00146       double incr=fval(ii)*xdiff/tmpsum;
00147       double fdiff = 1.-siglevel(level_index) - isum;
00148       if ( incr >= fdiff)
00149       {
00150         double delta=fdiff/incr;
00151         left_bd(level_index)=xs(ii)+delta*xdiff;
00152         break;
00153       }
00154       isum+=incr;
00155     }
00156     //cout << "tmpsum = " << tmpsum << endl;
00157     isum=0;
00158     for (ii=ub-1;ii>=lb;ii--)
00159     {
00160       xdiff=xs(ii+1)-xs(ii);
00161       double incr=fval(ii)*xdiff/tmpsum;
00162       double fdiff = 1.-siglevel(level_index) - isum;
00163       if ( incr >= fdiff)
00164       {
00165         double delta=fdiff/incr;
00166         right_bd(level_index)=xs(ii)+delta*xdiff;
00167         break;
00168       }
00169       isum+=incr;
00170     }
00171   }