clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% genera datos %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% model %%%%%%%%%%%%%%%%
%%%%% theta,t son variables genericas... de una funcion
f_gen=@(theta,t) sin(theta.*t);
%f_gen=@(theta,t) theta.*t;
%sig_r=1; %%% for simplicity we consider the std of the noise known
sig_r=0.5; %%% 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_gen(theta_true,t_datos)+sig_r*randn(1,Ndatos);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tc=0:0.001:10;
subplot(3,1,1)
plot(tc,f_gen(theta_true,tc),'k')
hold on
plot(t_datos,y,'ro')
xlabel('t')
set(gca,'Fontsize',15)
set(gca,'FontWeight','Bold')
title('data and model')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% since we are not "God" we could consider another model...
f=@(theta,t) sin(theta.*t);
%f=@(theta,t) theta.*t;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% forma mas rapida de calcular LIKE_tot 
%%LIKE_tot=@(theta) (1/(sqrt(2*pi*sig_r^2))).^Ndatos.*exp(-sum((y-f(theta)).^2)/(2*sig_r.^2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% theta,y,t son variables genericas... de una funcion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
like_piece=@(theta,y,t) 1/(sqrt(2*pi*sig_r^2)).*exp(-sum((y-f(theta,t)).^2)/(2*sig_r.^2));
step_ev=0.005;
theta_ev=-10:step_ev:10; %%% for plot
LIKE_tot=ones(1,length(theta_ev));
for j=1:length(theta_ev)  
          LIKE_tot(j)=LIKE_tot(j).*like_piece(theta_ev(j),y,t_datos);     
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,1,2)
plot(theta_ev,LIKE_tot)
title('Likelihood fuction')
xlabel('theta')
set(gca,'Fontsize',15)
xticks(gca,[-10 -7.5 -5 -4 -2 0 2 5 7.5 10])
set(gca,'FontWeight','Bold')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sig_p=2;
prior=@(theta) 1/sqrt(2*pi*sig_p^2)*exp(-theta.^2/(2*sig_p^2));
%par_unif=3;
%prior=@(theta) 1/(2*par_unif).*(theta>=-par_unif & theta<=par_unif);
POST=LIKE_tot.*prior(theta_ev);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,1,3)
plot(theta_ev,POST,'k')
title('Posterior')
xlabel('theta')
set(gca,'Fontsize',15)
xticks(gca,[-10 -7.5 -5 -4 -2 0 2 5 7.5 10])
set(gca,'FontWeight','Bold')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('----------------------------------------------------------------------')
disp('Riemman approximation of the Normalizing constant - marginal likelihood')
Zest=sum(POST*step_ev)
% % % disp('Riemman approximation of the first moment')
% % % mu_est=sum(theta_ev.*POST/Zest*step_ev)
% % % disp('Riemman approximation of the second moment')
% % % M2_est=sum(theta_ev.^2.*POST/Zest*step_ev)
% % % disp('Riemman approximation of the variance (central second moment)')
% % % VAR_est=sum((theta_ev-mu_est).^2.*POST/Zest*step_ev)
% % % %%% test %%%
% % % VAR_est=M2_est-(mu_est)^2
% % % disp('----------------------------------------------------------------------')
