% A matlab file for an exercise on Bayesian estimation
% Computation of MSE and LMSE estimates

clear
syms x y

%% Exercise 1.1
joint pdf
fxy=-3/2*x^2+2*x*y;
% estremi di integrazione
xl=0; xu=1;
yl=1; yu=2;

% %% Exercise 1.2
% fxy=(1/12)*(x+y)*exp(-y);
% % estremi di integrazione
% xl=0; xu=4;
% yl=0; yu=inf;

% marginal pdfs of x and y
fx=int(fxy,y,yl,yu);
fy=int(fxy,x,xl,xu);

% pdf of x conditioned to y
fxcy=fxy/fy;

% MSE estimate
xmse=int(x*fxcy,x,xl,xu)

% mean and covariances
mx=int(x*fx,x,xl,xu);
my=int(y*fy,y,yl,yu);
Rxy=int(int((x-mx)*(y-my)*fxy,y,yl,yu),x,xl,xu);
Ry=int((y-my)^2*fy,y,yl,yu);
Rx=int((x-mx)^2*fx,x,xl,xu);

% LMSE estimate
xlmse=mx+Rxy*inv(Ry)*(y-my)

% error variances
MSE=int(int((x-xmse)^2*fxy,y,yl,yu),x,xl,xu);
MSElinear=Rx-Rxy^2/Ry;
MSEprior=Rx;
% compare errors
double([MSE MSElinear MSEprior])

if yu==inf, yu=10; end
figure
ya=[yl:0.05:yu];
plot(ya,subs(xmse,ya),'r')
hold on
plot(ya,subs(xlmse,ya),'g')
plot(ya,mx*ones(size(ya)),'b')
xlabel('y')
ylabel('Estimates of X')
legend('MSE','LMSE','E[x]')