Gibbs Samplerと収束判定の例

#include <oxstd.h>
#include <oxprob.h>
#include <oxfloat.h>
#include <oxdraw.h>
//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));
}

//
main(){
decl a,b,Bm,burn,file,i,mu,mu_s,n,n_acf,n_x,sp,var,var_s,x,xbar;
// true values
 mu=1;var=2;
// generate x(1),...,x(n_x)
 n_x=100;x=mu+sqrt(var)*rann(n_x,1);xbar=meanc(x);
//Burn-in & Samples
 burn=1000;n=10000;
 mu_s=var_s=zeros(n,1);
 for(i=-burn;i<n;++i){
  mu=xbar+sqrt(var/n)*rann(1,1);
  // Generate var using Gamma(a,b)
  //  mean=a/b, variaince=a/b^2
  a=0.5*n_x;b=0.5*sumsqrc(x-mu);var=1.0/rangamma(1,1,a,b);
  if(i>=0){mu_s[i]=mu;var_s[i]=var;}
 }
//
println("XBAR=",xbar);
Bm=100;
println(
"%r",{"Mu"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke"},
meanc(mu_s)~varc(mu_s)^0.5~quantilec(mu_s,<0.025,0.975>)'~fGeweke(mu_s,Bm),"\n");
//
println("Std Err of Estimated Posterior Mean of Mu:");
println(" (1)Parzen Window: ",sqrt(fTsvar(mu_s,Bm)/n));
println(" (2)Batch Mean Method:",sqrt(fTsvar_batch(mu_s,Bm)/n));
//
println(
"%r",{"Sigma^2"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke"},
meanc(var_s)~varc(var_s)^0.5~quantilec(var_s,<0.025,0.975>)'~fGeweke(var_s,Bm),"\n");
//
 n_acf=50;
 //ACF
 DrawCorrelogram(0,mu_s',{"Mu"},n_acf);
 DrawTitle(0,"Sample ACF");
 DrawCorrelogram(1,var_s',{"Sigma^2"},n_acf);
 DrawTitle(1,"Sample ACF");
 //Path
 DrawTMatrix(2,mu_s',{"Mu"},1,1,1);
 DrawTitle(2,"Sample Path");
 DrawTMatrix(3,var_s',{"Sigma^2"},1,1,1);
 DrawTitle(3,"Sample Path");
 DrawDensity(4,mu_s',{"Mu"},1,0,0);
 //Posterior density
 DrawTitle(4,"Posterior Density");
 DrawDensity(5,var_s',{"Sigma^2"},1,1,0);
 DrawTitle(5,"Posterior Density");
 SaveDrawWindow("gibbs.ps");
 CloseDrawWindow();
// QQ-plot
decl t;
 t=(mu_s-xbar)/sqrt(varc(x)/n);
 DrawQQ(0,t',{"Mu"},QQ_T,n_x-1,0);
 DrawTitle(0,"QQ-plot vs T");
 DrawQQ(1,t',{"Mu"},QQ_N,0,0);
 DrawTitle(1,"QQ-plot vs Normal");
 SaveDrawWindow("qq.ps");
 CloseDrawWindow();

//
 file = fopen("mu.txt","w");fprint(file,"mu","%15.10f",mu_s);fclose(file);
 file = fopen("var.txt","w");fprint(file,"var","%15.10f",var_s);fclose(file);
 }