ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
newmodm3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: newmodm3.cpp 2472 2014-10-04 00:11:55Z 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(const int x, const int y)
00010 {
00011   return x < y ? y : x;
00012 }
00013 int xxxmin(const int x, const int y)
00014 {
00015   return x < y ? x : y;
00016 }
00020 void get_confidence_interval(const dvector& _left_bd, const dvector& _right_bd,
00021   dmatrix& ms, const dvector& xs, const dvector& siglevel,
00022   const int& level_index, int index)
00023   {
00024     dvector& left_bd=(dvector&) _left_bd;
00025     dvector& right_bd=(dvector&) _right_bd;
00026     dvector& fval=ms(index);
00027     int lb=0;
00028     int rb=0;
00029     double xdiff=0.0;
00030     double isum=0;
00031 
00032     int mmin=fval.indexmin();
00033     int mmax=fval.indexmax();
00034     double tmp=fval(mmin);
00035     int imax=mmin;
00036     int i;
00037     for (i=mmin+1;i<=mmax;i++)
00038     {
00039       if (fval(i)>tmp)
00040       {
00041         tmp=fval(i);
00042         imax=i;
00043       }
00044     }
00045     if (imax>mmin)
00046     {
00047       lb=imax-1;
00048       rb=imax;
00049     }
00050     else
00051     {
00052       lb=imax;
00053       rb=imax+1;
00054     }
00055     double integral=0.0;
00056 
00057     for (i=mmin+1;i<=mmax;i++)
00058     {
00059       integral+=fval(i)*(xs(i)-xs(i-1));
00060     }
00061     cout << integral << endl;
00062     for (;;)
00063     {
00064       if (lb <= fval.indexmin() &&  rb > fval.indexmax())
00065       {
00066         int lower=xxxmax(lb,fval.indexmin());
00067         int upper=xxxmin(rb-1,fval.indexmax());
00068         left_bd(level_index)=xs(lower);
00069         right_bd(level_index)=xs(upper);
00070         break;
00071       }
00072 
00073       if (fval(lb)>=fval(rb) || rb==fval.indexmax())
00074       {
00075         if (lb>fval.indexmin())
00076         {
00077           xdiff=xs(lb)-xs(lb-1);
00078         }
00079         else
00080         {
00081           int minind=fval.indexmin();
00082           xdiff=xs(minind+1)-xs(minind);
00083         }
00084         double incr=fval(lb)*xdiff/integral;
00085         //double incr=fval(lb);
00086         if ( incr >= (siglevel(level_index)-isum))
00087         {
00088           double diff = siglevel(level_index) - isum;
00089           double delta=diff/incr;
00090           if (lb>fval.indexmin())
00091           {
00092             left_bd(level_index)=xs(lb)-delta*(xdiff);
00093           }
00094           else
00095           {
00096             left_bd(level_index)=xs(fval.indexmin());
00097           }
00098           right_bd(level_index)=xs(rb);
00099           break;
00100         }
00101         isum+=incr;
00102         lb--;
00103       }
00104       else
00105       {
00106         xdiff=xs(rb)-xs(rb-1);
00107         double incr=fval(rb)*xdiff/integral;
00108         //double incr=fval(rb);
00109         if ( incr >= (siglevel(level_index)-isum) ||
00110           rb == mmax)
00111         {
00112           double diff = siglevel(level_index) - isum;
00113           double delta=diff/incr;
00114           left_bd(level_index)=xs(lb);
00115           right_bd(level_index)=xs(rb)+delta*(xdiff);
00116           break;
00117         }
00118         isum+=incr;
00119         rb++;
00120       }
00121     }
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   }