%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Metropolis Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Likelihood function or posterior distribution of normal mixture
%%% N(a,b^2) : Normal distribution, average a and variance b^2
%%% Model p=(1-aa) N(0,STD^2) + aa N(b,STD^2)
%%% True q=(1-AA)N(0,STD^2) + AA N(BB,STD^2)
%%% Uniform distribution is used as prior.
%%% You can rotate 3D figure by dragging the mouse.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;
clear;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NNN=100; %%% Sample Number %%% NNN=100, 200, ...
AA=0.2; %%% True mixture ratio
BB=1.5; %%% True b
RAND_1=1; %%% Random seed for mixture
RAND_2=1; %%% Random seed for normal distribution
STD=1.0; %%% standard deviation
BURNIN=2000; %%% Burn-in Number
MCMC=12000; %%% MCMC number
INTER=50; %%% MCMC interval
SIZEM=floor((MCMC-BURNIN)/INTER); %%% MCMC sample number
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sample generate %%%%%%%
rand('state',RAND_1);
randn('state',RAND_2);
XX=STD*randn(1,NNN);
YY=rand(1,NNN);
for i=1:1:NNN
if(YY(i)1)
aa=aa-1;
end
DD=HHH(aa,bb)-HHH(aa0,bb0);
if(exp(-DD)>rand(1,1))
aa0=aa;
bb0=bb;
change=change+1;
end
if(mod(t,INTER)==1&t>BURNIN)
k=k+1;
waa(k)=aa0;
wbb(k)=bb0;
MCN(k)=t;
EEE(k)=HHH(aa0,bb0);
end
end
%%%%%%%%MCMC end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MCMC_exchange_rate=change/MCMC
subplot(1,2,1);
plot(waa,wbb,'r+'); hold on
xlabel('Mixture Ratio');
ylabel('parameter b');
subplot(1,2,2);
plot(MCN,EEE,'b');
xlabel('MCMC TIMES');
ylabel('Energy'); hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%