stats.m.txt

这是马尔可夫-蒙特卡罗算法的MATLAB源程序.

function [fbar,stdf,tauf]=stats(v,mx,StatNames)

[s,N]=size(v);

for f=1:s
   w=v(f,1:N);
   r=zeros(1,mx+1); 
   N=length(w);
   for k=0:mx
      CC=corrcoef([w(1:(N-k))',w((1+k):N)']);
      r(k+1)=CC(1,2);
   end

   G=zeros(1,mx);
   for k=1:mx
      G(k)=r(k+1)+r(k);
   end
   
   tauf(f)=-r(1)+2*G(1);
   avt=1;
   for M=1:(mx-1)
      if G(M+1)<G(M) & G(M+1)>0
         tauf(f)=tauf(f)+2*G(M+1);
         avt=avt+2*r(M+1)^2;
      else
         break;
      end
   end
   if tauf(f)<1, tauf(f)=1; if M~=1, disp('unexpected M>1 in stats');keyboard;pause; end, end %when r(2)<0
   avt=avt/N;
   ast=sqrt(avt);
   
   fbar(f)=mean(w);
   stdf(f)=sqrt(tauf(f)/N)*std(w);
   
%   subplot(s,2,2*f-1);
   subplot(2,2,2*f-1);
   plot(0:mx,r,[0,mx],[0,0],'r',[0,mx],[2*ast,2*ast],'r--',[0,mx],[-2*ast,-2*ast],'r--');
   if nargin==3
       %      title([StatNames{f},' \tau_f=',num2str(tauf(f)),' M=',num2str(M),' N=',num2str(N)]);
   else
      title([' \tau_f=',num2str(tauf(f)),' M=',num2str(M),' N=',num2str(N)]);
   end
   ylabel('autocorrelation');
%   subplot(s,2,2*f);
  subplot(2,2,2*f);
   plot(w);
   xy=axis;xy(2)=N;
   axis(xy);
   if nargin==3
      ylabel(StatNames{f});
   end
end
%subplot(s,2,2);
subplot(2,2,2);
%title('MCMC output trace');
%subplot(s,2,2*s-1);
subplot(2,2,2*s-1);
xlabel('lag');
%subplot(s,2,2*s);
subplot(2,2,2*s);
xlabel('MCMC updates');