ADMB Documentation  11.1.2438
 All Classes Files Functions Variables Typedefs Friends Defines
c_ghk.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: c_ghk.cpp 1678 2014-02-25 21:31:02Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00017 double ghk(const dvector& lower,const dvector& upper,const dmatrix& Sigma,
00018   const dmatrix& eps)
00019 {
00020   int m=eps.indexmax();
00021   int n=lower.indexmax();
00022   double ssum=0.0;
00023   dmatrix ch=choleski_decomp(Sigma);
00024   dvector l(1,n);
00025   dvector u(1,n);
00026 
00027   for (int k=1;k<=m;k++)
00028   {
00029     double weight=1.0;
00030     l=lower;
00031     u=upper;
00032     for (int j=1;j<=n;j++)
00033     {
00034       l(j)/=ch(j,j);
00035       u(j)/=ch(j,j);
00036       double Phiu=cumd_norm(u(j));
00037       double Phil=cumd_norm(l(j));
00038       weight*=Phiu-Phil;
00039       double eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil);
00040       for (int i=j+1;i<=n;i++)
00041       {
00042         double tmp=ch(i,j)*eta;
00043         l(i)-=tmp;
00044         u(i)-=tmp;
00045       }
00046     }
00047     ssum+=weight;
00048   }
00049   return ssum/m;
00050 }
00051 
00056 void ghk_test(const dmatrix& eps,int i)
00057 {
00058   if (i<eps.indexmin())
00059   {
00060     cerr << "Index too low in function ghk -- min is "
00061          << eps.indexmin() << " you have " << i << endl;
00062     exit(21);
00063   }
00064   else  if (i>eps.indexmax())
00065   {
00066     cerr << "Index too high in function ghk -- max is "
00067          << eps.indexmax() << " you have " << i << endl;
00068     exit(21);
00069   }
00070 }
00071 
00076 double ghk(const dvector& lower,const dvector& upper,const dmatrix& Sigma,
00077   const dmatrix& eps,int i)
00078 {
00079   int n=lower.indexmax();
00080   dmatrix ch=choleski_decomp(Sigma);
00081   dvector l(1,n);
00082   dvector u(1,n);
00083 
00084   ghk_test(eps,i); // test for valid i range
00085   double weight=1.0;
00086   int k=i;
00087   {
00088     l=lower;
00089     u=upper;
00090     for (int j=1;j<=n;j++)
00091     {
00092       l(j)/=ch(j,j);
00093       u(j)/=ch(j,j);
00094       double Phiu=cumd_norm(u(j));
00095       double Phil=cumd_norm(l(j));
00096       weight*=Phiu-Phil;
00097       double eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil);
00098       for (int i=j+1;i<=n;i++)
00099       {
00100         double tmp=ch(i,j)*eta;
00101         l(i)-=tmp;
00102         u(i)-=tmp;
00103       }
00104     }
00105   }
00106   return weight;
00107 }