clear all
close all
clc
%%%%%%%%%%%%%%%
y=5; %%% data 
s=2; %%% std noise
%%%%%%%%%%%%%%%
G=@(theta) theta.^2; %%% OBervation MODEL 
%G=@(theta) sin(theta);
%%%% POST FOR THE MH ALGORITHM REPRESENTS ANY TARGET DENSITIES...
Post=@(theta) exp(-(y-G(theta)).^2/(2*s.^2));
s_prop=3;
Prop=@(v1,v2) exp(-(v1-v2).^2/(2*s_prop^2));
%%%%%%%%
theta=-10:0.01:10;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(theta,Post(theta),'b')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta_c=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=1000;
%%%%%%%%%%%%%%%%%
for t=2:T+1
    theta_pr=theta_c(t-1)+s_prop*randn(1,1); %%% proposal q(theta|theta_{t-1})
    Tar_now=Post(theta_pr);
    Tar_bef=Post(theta_c(t-1));
    Prop_now=Prop(theta_pr,theta_c(t-1));
    Prop_bef=Prop(theta_c(t-1),theta_pr);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    rho=(Tar_now*Prop_bef)/(Tar_bef*Prop_now);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    a=min([1,rho]); %%% acceptance function
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    u=rand(1,1);
    %%%%%%%%%%%%%%
    if u<=a
        theta_c(t)=theta_pr;
    else
        theta_c(t)=theta_c(t-1);
    end
  %%%%%%%  
end
figure
hist(theta_c,30)
%%%%%%%%%%%%%%%%%%%%
Zapp=sum(Post(theta).*0.01)
mu_est=sum(theta.*1/Zapp.*Post(theta).*0.01)
var_est=sum((theta-mu_est).^2.*1/Zapp.*Post(theta).*0.01)
var(theta_c)
%%%%%%%%%%%%%%%%%%%%
figure
plot(theta,1/Zapp.*Post(theta),'k','LineWidth',5)
%%%%%%%%%%%%%%%%%%%%
[d1,d2]=hist(theta_c,30);
step_now=d2(2)-d2(1);
ZappHist=sum(d1.*step_now);
hold on
bar(d2,1/ZappHist*d1,'r')
plot(theta,1/Zapp.*Post(theta),'k','LineWidth',5)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(theta_c(1:100))


