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')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu_prop=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=200000;
%%%%%%%%%%%%%%%%%
  theta_pr=mu_prop+s_prop*randn(1,N); %%% 
  for j=1:length(theta_pr)  
          Tar_now(j)=Post(theta_pr(j));     
  end
   Prop_now=Prop(theta_pr,mu_prop);
   w=Tar_now./Prop_now;
   wn=w./sum(w);
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%% end %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('estimation with IS')
Z_est_IS= mean(w)
mu_est_IS= sum(wn.*theta_pr)
M2_est_IS= sum(wn.*theta_pr.^2)
var_est_IS= M2_est_IS-(mu_est_IS).^2
disp('----------------------------------------------------------------------')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



