MHアルゴリズムの例
プログラム
#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 burn,i,n,nobs,phi,phi_n,phi_o,phi_s,mu_phi,var_phi,var_e,y; decl Bm,file,n_acf; nobs=1001;y=zeros(nobs,1);phi=0.9;var_e=1.0; // Generate Dataset to analyze y[0]=rann(1,1)*sqrt(var_e/(1-phi^2)); for(i=0;i<nobs-1;++i){ y[i+1]=phi*y[i]+rann(1,1)*sqrt(var_e); } // Dataset // Prior of phi ~ Uniform (-1,1), f(phi)=0.5 on (-1,1) n=10000;burn=0.1*n;phi_s=zeros(n,1); for(i=-burn;i<n;++i){ phi_o=phi; mu_phi=(y[0:nobs-2]'*y[1:nobs-1])/sumsqrc(y[1:nobs-2]); var_phi=var_e/sumsqrc(y[1:nobs-2]); // Generate phi ~ TN(mu_phi,var_phi) on (-1,1) do{ phi_n=mu_phi+rann(1,1)*sqrt(var_phi); }while(fabs(phi_n)>=1); // MH if(ranu(1,1)<= (0.5*sqrt(1-phi_n^2))/(0.5*sqrt(1-phi_o^2))){phi=phi_n;} else{phi=phi_o;} if(i>=0){phi_s[i]=phi;} } // true values Bm=0.1*n; println( "%r",{"PHI"}, "%c",{"Mean","Stdev","95%L","95%U","Geweke"}, meanc(phi_s)~varc(phi_s)^0.5~quantilec(phi_s,<0.025,0.975>)'~fGeweke(phi_s,Bm),"\n"); // n_acf=50; //ACF DrawCorrelogram(0,phi_s',{"PHI"},n_acf); DrawTitle(0,"Sample ACF"); //Path DrawTMatrix(1,phi_s',{"PHI"},1,1,1); DrawTitle(1,"Sample Path"); //Posterior density DrawDensity(2,phi_s',{"PHI"},1,0,0); DrawTitle(2,"Posterior Density"); SaveDrawWindow("mh.eps"); CloseDrawWindow(); // file = fopen("phi.txt","w");fprint(file,"phi","%15.10f",phi_s);fclose(file); }
結果
Mean Stdev 95%L 95%U Geweke PHI 0.89917 0.014046 0.87210 0.92694 0.99415