df1b2im5.cpp
Go to the documentation of this file.
```00001 /*
00002  * \$Id: df1b2im5.cpp 1136 2013-08-05 19:40:39Z johnoel \$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #if defined(USE_LAPLACE)
00013 #  include <df1b2fun.h>
00015
00016 static void xxx(void){;}
00017
00022 double calculate_importance_sample_block_diagonal_option_antithetical
00023   (const dvector& x,const dvector& u0,const dmatrix& Hess,
00025   const dmatrix& _Hessadjoint,function_minimizer * pmin)
00026 {
00030   const int xs=x.size();
00031   const int us=u0.size();
00033   int nsc=pmin->lapprox->num_separable_calls;
00034   const ivector lrea = (*pmin->lapprox->num_local_re_array)(1,nsc);
00035   int hroom =  sum(square(lrea));
00036   int nvar=x.size()+u0.size()+hroom;
00037   independent_variables y(1,nvar);
00038
00039   // need to set random effects active together with whatever
00040   // init parameters should be active in this phase
00041   initial_params::set_inactive_only_random_effects();
00042   initial_params::set_active_random_effects();
00043   /*int onvar=*/initial_params::nvarcalc();
00044   initial_params::xinit(y);    // get the initial values into the
00045   // do we need this next line?
00046   y(1,xs)=x;
00047
00048   int i,j;
00049
00050   // contribution for quadratic prior
00052   {
00054     int & vxs = (int&)(xs);
00056   }
00057  // Here need hooks for sparse matrix structures
00058
00059   dvar3_array & block_diagonal_vhessian=
00060     *pmin->lapprox->block_diagonal_vhessian;
00061   block_diagonal_vhessian.initialize();
00062   dvar3_array& block_diagonal_ch=
00063     *pmin->lapprox->block_diagonal_vch;
00064     //dvar3_array(*pmin->lapprox->block_diagonal_ch);
00065   int ii=xs+us+1;
00066   d3_array& bdH=(*pmin->lapprox->block_diagonal_hessian);
00067   int ic;
00068   for (ic=1;ic<=nsc;ic++)
00069   {
00070     int lus=lrea(ic);
00071     for (i=1;i<=lus;i++)
00072       for (j=1;j<=lus;j++)
00073         y(ii++)=bdH(ic)(i,j);
00074   }
00075
00076   dvector g(1,nvar);
00079   dvar_vector vy=dvar_vector(y);
00080   //initial_params::stddev_vscale(d,vy);
00081   ii=xs+us+1;
00082   if (initial_df1b2params::have_bounded_random_effects)
00083   {
00084     cerr << "can't do importance sampling with bounded random effects"
00085      " at present" << endl;
00087   }
00088   else
00089   {
00090     int ic;
00091     for (ic=1;ic<=nsc;ic++)
00092     {
00093       int lus=lrea(ic);
00094       for (i=1;i<=lus;i++)
00095       {
00096         for (j=1;j<=lus;j++)
00097         {
00098           block_diagonal_vhessian(ic,i,j)=vy(ii++);
00099         }
00100       }
00101       block_diagonal_ch(ic)=
00102         choleski_decomp(inv(block_diagonal_vhessian(ic)));
00103     }
00104   }
00105
00106    int nsamp=pmin->lapprox->num_importance_samples;
00107
00108    dvariable vf=0.0;
00109
00110    dvar_vector sample_value(1,nsamp);
00111    sample_value.initialize();
00112
00113    dvar_vector tau(1,us);;
00114    int is;
00115    for (is=1;is<=nsamp;is++)
00116    {
00117      int offset=0;
00118      pmin->lapprox->importance_sampling_counter=is;
00119      int ic;
00120      for (ic=1;ic<=nsc;ic++)
00121      {
00122        int lus=lrea(ic);
00123        tau(offset+1,offset+lus).shift(1)=block_diagonal_ch(ic)*
00124          (*pmin->lapprox->antiepsilon)(is);
00125        offset+=lus;
00126      }
00127
00128      // have to reorder the terms to match the block diagonal hessian
00129      imatrix & ls=*(pmin->lapprox->block_diagonal_re_list);
00130      int mmin=ls.indexmin();
00131      int mmax=ls.indexmax();
00132
00133      int ii=1;
00134      int i;
00135      for (i=mmin;i<=mmax;i++)
00136      {
00137        int cmin=ls(i).indexmin();
00138        int cmax=ls(i).indexmax();
00139        for (int j=cmin;j<=cmax;j++)
00140        {
00141          vy(ls(i,j))+=tau(ii++);
00142        }
00143      }
00144      if (ii-1 != us)
00145      {
00146        cerr << "error in interface" << endl;
00148      }
00149      initial_params::reset(vy);    // get the values into the model
00150      ii=1;
00151      for (i=mmin;i<=mmax;i++)
00152      {
00153        int cmin=ls(i).indexmin();
00154        int cmax=ls(i).indexmax();
00155        for (int j=cmin;j<=cmax;j++)
00156        {
00157          vy(ls(i,j))-=tau(ii++);
00158        }
00159      }
00160
00161      *objective_function_value::pobjfun=0.0;
00162      //int istop=0;
00163      if (is==65)
00164      {
00165         xxx();
00166      }
00168
00169      if (pmin->lapprox->use_outliers==0)
00170      {
00171        // assumes that all separable calls have the same number
00172        // of random effects
00173        double neps=0.5*nsc*norm2((*pmin->lapprox->antiepsilon)(is));
00174
00175        (*pmin->lapprox->importance_sampling_values)(is)=
00176            neps-value(*objective_function_value::pobjfun);
00177
00178        (*pmin->lapprox->importance_sampling_weights)(is)=neps;
00179
00180         sample_value(is)=*objective_function_value::pobjfun
00181           -neps;
00182
00183      }
00184      else
00185      {
00186        dvector& e=pmin->lapprox->epsilon(is);
00187        double neps=-sum(log(.95*exp(-0.5*square(e))+
00188           0.0166666667*exp(-square(e)/18.0)));
00189
00190        (*pmin->lapprox->importance_sampling_values)(is)=
00191          neps-value(*objective_function_value::pobjfun);
00192
00193        sample_value(is)=*objective_function_value::pobjfun
00194          -neps;
00195      }
00196    }
00197
00198    nsc=pmin->lapprox->num_separable_calls;
00199    dmatrix weights(1,nsc,1,nsamp);
00200    for (is=1;is<=nsamp;is++)
00201    {
00202      int offset=0;
00203      int ic;
00204      for (ic=1;ic<=nsc;ic++)
00205      {
00206        int lus=lrea(ic);
00207        // assumes that all spearable calls have the same number of
00208        // random effects
00209        dvector e= (*pmin->lapprox->antiepsilon)(is);
00210        offset+=lus;
00211        if (pmin->lapprox->use_outliers==0)
00212        {
00213          weights(ic,is)=0.5*norm2(e);
00214        }
00215        else
00216        {
00217          weights(ic,is)=-sum(log(.95*exp(-0.5*square(e))+
00218           0.0166666667*exp(-square(e)/18.0)));
00219        }
00220      }
00221    }
00222    dvar_vector lcomp(1,nsc);
00223    for (int isc=1;isc<=nsc;isc++)
00224    {
00225      dvar_vector & comps=
00226        (*pmin->lapprox->importance_sampling_components)(isc);
00227
00228      dvar_vector diff=comps-weights(isc);
00229      double dmin=min(value(diff));
00230      lcomp(isc)=dmin-log(mean(exp(dmin-diff)));
00231    }
00232
00233    //double ns=lcomp.indexmax()-lcomp.indexmin()+1;
00234    //double min_vf=min(value(lcomp));
00235    vf= sum(lcomp);
00236    vf-=us*0.91893853320467241;
00237
00238    int sgn=0;
00239    dvariable ld=0.0;
00241    {
00242      for (int ic=1;ic<=nsc;ic++)
00243      {
00244        ld+=ln_det(block_diagonal_vhessian(ic),sgn);
00245      }
00246      ld*=0.5;
00247    }
00248    else
00249    {
00250      for (int ic=1;ic<=nsc;ic++)
00251      {
00252        ld+=ln_det_choleski(block_diagonal_vhessian(ic));
00253      }
00254      ld*=0.5;
00255    }
00256
00257    vf+=ld;
00258
00259    double f=value(vf);
00261
00262    // put uhat back into the model
00264    vy(xs+1,xs+us).shift(1)=u0;
00265    initial_params::reset(vy);    // get the values into the model
00267
00268   ii=1;
00269   for (i=1;i<=xs;i++)
00271   for (i=1;i<=us;i++)
00273   for (ic=1;ic<=nsc;ic++)
00274   {
00275     int lus=lrea(ic);
00276     for (i=1;i<=lus;i++)
00277     {
00278       for (j=1;j<=lus;j++)
00279       {