clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% genera datos %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% model %%%%%%%%%%%%%%%%
%%%%% theta,t son variables genericas... de una funcion
f=@(theta,t) sin(theta.*t);
%f=@(theta,t) theta.*t;
%sig_r=1; %%% for simplicity we consider the std of the noise known
sig_r=1; %%% for simplicity we consider the std of the noise known
%%%% datos observado %%%%
theta_true=2;
t_datos=0:1:10;
Ndatos=length(t_datos);
y=f(theta_true,t_datos)+sig_r*randn(1,Ndatos);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sig_p=2;
theta_ev=2;
Post=@(theta) Evaluate_Post(theta,y,t_datos,f,sig_r,sig_p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% compute groundtruth %%%%
step_ev=0.01;
theta_ev=-10:step_ev:10;
for j=1:length(theta_ev)  
          POST_tot_ev(j)=Post(theta_ev(j));     
end
disp('Riemman approximation of the Normalizing constant - marginal likelihood')
Zest=sum(POST_tot_ev*step_ev)
 disp('Riemman approximation of the first moment')
 mu_est=sum(theta_ev.*POST_tot_ev/Zest*step_ev)
 disp('Riemman approximation of the second moment')
 M2_est=sum(theta_ev.^2.*POST_tot_ev/Zest*step_ev)
disp('Riemman approximation of the variance (central second moment)')
 VAR_est=sum((theta_ev-mu_est).^2.*POST_tot_ev/Zest*step_ev)
 %%% test %%%
%%% VAR_est=M2_est-(mu_est)^2
 disp('----------------------------------------------------------------------')
%return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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=100000;
%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%% end MH %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% plot post %%%%%
step_ev=0.005;
theta_ev=-10:step_ev:10; %%% for plot
for j=1:length(theta_ev)  
          Post_tot(j)=Post(theta_ev(j));     
end
Zest=sum(Post_tot.*step_ev);
figure
[d1,d2]=hist(theta_c,70);
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)
plot(theta_ev,1/Zest.*Post_tot,'k','LineWidth',2)
grid on
figure
 plot(theta_c(1:100))
%%%%%%%%%%%%%%%%%%%%
%%%%
disp('estimation with MH')
mu_est_mh= mean(theta_c)
M2_est_mh= mean(theta_c.^2)
var_est_mh= var(theta_c)
disp('----------------------------------------------------------------------')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



