ADMB Documentation  11.1.2497
 All Classes Files Functions Variables Typedefs Friends Defines
v_ghk.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: v_ghk.cpp 1708 2014-02-28 21:43:41Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00017 dvariable ghk(const dvar_vector& lower,const dvar_vector& upper,
00018   const dvar_matrix& Sigma, const dmatrix& eps)
00019 {
00020   RETURN_ARRAYS_INCREMENT();
00021   int n=lower.indexmax();
00022   int m=eps.indexmax();
00023   dvariable ssum=0.0;
00024   dvar_matrix ch=choleski_decomp(Sigma);
00025   dvar_vector l(1,n);
00026   dvar_vector u(1,n);
00027 
00028   for (int k=1;k<=m;k++)
00029   {
00030     dvariable weight=1.0;
00031     l=lower;
00032     u=upper;
00033     for (int j=1;j<=n;j++)
00034     {
00035       l(j)/=ch(j,j);
00036       u(j)/=ch(j,j);
00037       dvariable Phiu=cumd_norm(u(j));
00038       dvariable Phil=cumd_norm(l(j));
00039       weight*=Phiu-Phil;
00040       dvariable eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil+1.e-30);
00041       for (int i=j+1;i<=n;i++)
00042       {
00043         dvariable tmp=ch(i,j)*eta;
00044         l(i)-=tmp;
00045         u(i)-=tmp;
00046       }
00047     }
00048     ssum+=weight;
00049   }
00050   RETURN_ARRAYS_DECREMENT();
00051   return ssum/m;
00052 }
00053 
00058 dvariable ghk_m(const dvar_vector& upper,const dvar_matrix& Sigma,
00059   const dmatrix& eps)
00060 {
00061   RETURN_ARRAYS_INCREMENT();
00062   int n=upper.indexmax();
00063   int m=eps.indexmax();
00064   dvariable ssum=0.0;
00065   dvar_vector u(1,n);
00066   dvar_matrix ch=choleski_decomp(Sigma);
00067 
00068   for (int k=1;k<=m;k++)
00069   {
00070     dvariable weight=1.0;
00071     u=upper;
00072     for (int j=1;j<=n;j++)
00073     {
00074       u(j)/=ch(j,j);
00075       dvariable Phiu=cumd_norm(u(j));
00076       weight*=Phiu;
00077       dvariable eta=inv_cumd_norm(1.e-30+Phiu*eps(k,j));
00078       for (int i=j+1;i<=n;i++)
00079       {
00080         dvariable tmp=ch(i,j)*eta;
00081         u(i)-=tmp;
00082       }
00083     }
00084     ssum+=weight;
00085   }
00086   RETURN_ARRAYS_DECREMENT();
00087   return ssum/m;
00088 }
00089 
00094 dvariable ghk_choleski(const dvar_vector& lower,const dvar_vector& upper,
00095   const dvar_matrix& ch, const dmatrix& eps)
00096 {
00097   RETURN_ARRAYS_INCREMENT();
00098   int n=lower.indexmax();
00099   int m=eps.indexmax();
00100   dvariable ssum=0.0;
00101   dvar_vector l(1,n);
00102   dvar_vector u(1,n);
00103 
00104   for (int k=1;k<=m;k++)
00105   {
00106     dvariable weight=1.0;
00107     l=lower;
00108     u=upper;
00109     for (int j=1;j<=n;j++)
00110     {
00111       l(j)/=ch(j,j);
00112       u(j)/=ch(j,j);
00113       dvariable Phiu=cumd_norm(u(j));
00114       dvariable Phil=cumd_norm(l(j));
00115       weight*=Phiu-Phil;
00116       dvariable eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil);
00117       for (int i=j+1;i<=n;i++)
00118       {
00119         dvariable tmp=ch(i,j)*eta;
00120         l(i)-=tmp;
00121         u(i)-=tmp;
00122       }
00123     }
00124     ssum+=weight;
00125   }
00126   RETURN_ARRAYS_DECREMENT();
00127   return ssum/m;
00128 }
00129 
00134 dvariable ghk_choleski_m(const dvar_vector& upper,
00135   const dvar_matrix& ch, const dmatrix& eps)
00136 {
00137   RETURN_ARRAYS_INCREMENT();
00138   int n=upper.indexmax();
00139   int m=eps.indexmax();
00140   dvariable ssum=0.0;
00141   dvar_vector u(1,n);
00142 
00143   for (int k=1;k<=m;k++)
00144   {
00145     dvariable weight=1.0;
00146     u=upper;
00147     for (int j=1;j<=n;j++)
00148     {
00149       u(j)/=ch(j,j);
00150       dvariable Phiu=cumd_norm(u(j));
00151       weight*=Phiu;
00152       dvariable eta=inv_cumd_norm(1.e-30+Phiu*eps(k,j));
00153       for (int i=j+1;i<=n;i++)
00154       {
00155         dvariable tmp=ch(i,j)*eta;
00156         u(i)-=tmp;
00157       }
00158     }
00159     ssum+=weight;
00160   }
00161   RETURN_ARRAYS_DECREMENT();
00162   return ssum/m;
00163 }
00164 
00165 void ghk_test(const dmatrix& eps,int i);
00166 
00171 dvariable ghk(const dvar_vector& lower,const dvar_vector& upper,
00172   const dvar_matrix& Sigma, const dmatrix& eps,int i)
00173 {
00174   RETURN_ARRAYS_INCREMENT();
00175   int n=lower.indexmax();
00176   dvar_matrix ch=choleski_decomp(Sigma);
00177   dvar_vector l(1,n);
00178   dvar_vector u(1,n);
00179 
00180   ghk_test(eps,i);
00181 
00182   dvariable weight=1.0;
00183   int k=i;
00184   {
00185     l=lower;
00186     u=upper;
00187     for (int j=1;j<=n;j++)
00188     {
00189       l(j)/=ch(j,j);
00190       u(j)/=ch(j,j);
00191       dvariable Phiu=cumd_norm(u(j));
00192       dvariable Phil=cumd_norm(l(j));
00193       weight*=Phiu-Phil;
00194       dvariable eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil+1.e-30);
00195       for (int i=j+1;i<=n;i++)
00196       {
00197         dvariable tmp=ch(i,j)*eta;
00198         l(i)-=tmp;
00199         u(i)-=tmp;
00200       }
00201     }
00202   }
00203   RETURN_ARRAYS_DECREMENT();
00204   return weight;
00205 }
00206 
00211 dvariable ghk_choleski_m_cauchy(const dvar_vector& upper,
00212   const dvar_matrix& ch, const dmatrix& eps)
00213 {
00214   RETURN_ARRAYS_INCREMENT();
00215   int n=upper.indexmax();
00216   int m=eps.indexmax();
00217   dvariable ssum=0.0;
00218   dvar_vector u(1,n);
00219 
00220   for (int k=1;k<=m;k++)
00221   {
00222     dvariable weight=1.0;
00223     u=upper;
00224     for (int j=1;j<=n;j++)
00225     {
00226       u(j)/=ch(j,j);
00227       dvariable Phiu=cumd_cauchy(u(j));
00228       weight*=Phiu;
00229       dvariable eta=inv_cumd_cauchy(Phiu*eps(k,j));
00230       for (int i=j+1;i<=n;i++)
00231       {
00232         dvariable tmp=ch(i,j)*eta;
00233         u(i)-=tmp;
00234       }
00235     }
00236     ssum+=weight;
00237   }
00238   RETURN_ARRAYS_DECREMENT();
00239   return ssum/m;
00240 }
00241 
00246 dvariable ghk_choleski_m_logistic(const dvar_vector& upper,
00247   const dvar_matrix& ch, const dmatrix& eps)
00248 {
00249   RETURN_ARRAYS_INCREMENT();
00250   int n=upper.indexmax();
00251   int m=eps.indexmax();
00252   dvariable ssum=0.0;
00253   dvar_vector u(1,n);
00254 
00255   for (int k=1;k<=m;k++)
00256   {
00257     dvariable weight=1.0;
00258     u=upper;
00259     for (int j=1;j<=n;j++)
00260     {
00261       u(j)/=ch(j,j);
00262       dvariable Phiu=cumd_logistic(u(j));
00263       weight*=Phiu;
00264       dvariable eta=inv_cumd_logistic(Phiu*eps(k,j));
00265       for (int i=j+1;i<=n;i++)
00266       {
00267         dvariable tmp=ch(i,j)*eta;
00268         u(i)-=tmp;
00269       }
00270     }
00271     ssum+=weight;
00272   }
00273   RETURN_ARRAYS_DECREMENT();
00274   return ssum/m;
00275 }