AR-MHアルゴリズム

プログラム

#include <oxstd.h>
#include <oxprob.h>
#include <oxfloat.h>
#include <oxdraw.h>
#import<maximize>
//Calculate variance of time series 
fTsvar(const x, const Bm){
decl n,sp;
// Calculate Spectral density function using
// Bandwidth=Bm at 2 points (0, pi) and Parzen  Window
n=rows(x);sp=periodogram(x,Bm,2,1)/n;
return M_2PI*sp[0];
}
fTsvar_batch(const x, const n){
// N=n*m
decl i,m,N,y,xbar;
N=rows(x);m=N/n;xbar=meanc(x);y=zeros(m,1);
for(i=0;i<m;++i){y[i]=meanc(x[i*n:(i+1)*n-1]);}
return n/(m-1)*sumsqrc(y-xbar);
}
// Calculate p-value for Convergence (Geweke's Method)
// H0: Convergence, H1: No Convergence
fGeweke(const x, const Bm){
decl n,n1,n2,sp1,sp2,var1,var2,x1,x2,x1bar,x2bar,z;
 n=rows(x);n1=floor(0.1*n);n2=floor(0.5*n);
 x1=x[:n1-1];x2=x[n-n2:];
 x1bar=meanc(x1);x2bar=meanc(x2);
 var1=fTsvar(x1,Bm);var2=fTsvar(x2,Bm);
 z=(x1bar-x2bar)/sqrt(var1/n1+var2/n2);
return 2*tailn(fabs(z));
}


static decl x;
fLk(const dA, const adFunc, const avScore, const amHess)
{
decl alpha,nobs;nobs=rows(x);alpha=dA;
adFunc[0]=-nobs*loggamma(alpha)+(alpha-1)*sumc(log(x))-sumc(x)-log(alpha);
if(avScore){avScore[0]=-nobs*polygamma(alpha,0)+sumc(log(x))-1/alpha;}
if(amHess){amHess[0]=-nobs*polygamma(alpha,1)+1/alpha^2;}
return 1;
}
//
main(){
decl a,a_m,a_n,a_o,a_s,i,n,nobs;
decl frac,h_n,h_o,l_m,l_n,l_o,l1_m,l2_m,mu,var;
decl burn,ar_count,ar_accept,mh_accept,mh_count;
decl Bm,file,n_acf;
a=1;nobs=1000;
// Generate Dataset X_1,...,X_n ~ Gamma(a,1)
x=rangamma(nobs,1,a,1);
// Find the mode of L(a). initial value: a
  a_m=a;MaxBFGS(fLk,&a_m,&l_m,0,0);println("mode=",a_m);
  fLk(a_m,&l_m,&l1_m,&l2_m);var=-1/l2_m;mu=a_m+l1_m*var;
//
n=10000;burn=0.1*n;a_s=zeros(n,1);
ar_count=ar_accept=mh_count=mh_accept=0;
for(i=-burn;i<n;++i){
   a_o=a;fLk(a_o,&l_o,0,0);
   // AR Step
   if(fmod(i,100)==0){println(i~a);}
   h_o=l_m+(a_o-a_m)*l1_m+0.5*(a_o-a_m)^2*l2_m;
   do{
    do{a_n=mu+rann(1,1)*sqrt(var);}while(a_n<=0);  
        fLk(a_n,&l_n,0,0);
    h_n=l_m+(a_n-a_m)*l1_m+0.5*(a_n-a_m)^2*l2_m;
        ar_count+=1;
   } while(ranu(1,1)>=exp(l_n-h_n));
    ar_accept+=1;
 // MH Step     
  frac=exp(l_n+min(l_o|h_o)-l_o-min(l_n|h_n));
  if(ranu(1,1) <= frac){a=a_n;mh_accept+=1;}
  mh_count+=1;
  if(i>=0){a_s[i]=a;}
}
println("AR acceptance=",ar_accept/ar_count);
println("MH acceptance=",mh_accept/mh_count);
Bm=0.1*n;
//Bm=0;
println(
"%r",{"ALPHA"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke"},
meanc(a_s)~varc(a_s)^0.5~quantilec(a_s,<0.025,0.975>)'~fGeweke(a_s,Bm),"\n");
//
 n_acf=500;
 //ACF
 DrawCorrelogram(0,a_s',{"ALPHA"},n_acf);
 DrawTitle(0,"Sample ACF");
 //Path
 DrawTMatrix(1,a_s',{"ALPHA"},1,1,1);
 DrawTitle(1,"Sample Path");
 //Posterior density
 DrawDensity(2,a_s',{"ALPHA"},1,1,0);
 DrawTitle(2,"Posterior Density");
 SaveDrawWindow("ar-mh.ps");
 CloseDrawWindow();
//
 file = fopen("a.txt","w");fprint(file,"a","%15.10f",a_s);fclose(file);
}

結果

mode=
       1.0090

      -1000.0       1.0000

      -900.00       1.0080
　　　（途中略）

       9600.0       1.0030

       9700.0       1.0262

       9800.0       1.0400

       9900.0      0.98925
AR acceptance=0.99719
MH acceptance=0.995818

               Mean        Stdev         95%L         95%U       Geweke
ALPHA        1.0092     0.024840      0.96075       1.0586     0.090136