% Add these directory to the search path of MATLAB
clear; clear functions; clf;

disp(' PROGRAM to compare Least Square estimate with the ones obtained');
disp(' using the Joint PDF.');
disp(' ');

% Generate two random variable Y and X with a non-zero linear relationship
% Choose example number from class.
example=input('Choose example Number (there are 4 examples total)   => ');
% Choose number of samples in the generation of random variables.
NS=10000;
% Call function to generate the variables.
[x,y,n]=bs_generate_XY(NS, example);

% Keep NS_indep number of independent sample for model
% testing.
NS_indep=50;
x_indep=x(1:NS_indep);
y_indep=y(1:NS_indep);

% Take this sub-sample out of the X and Y data
x=x(NS_indep+1:end);
y=y(NS_indep+1:end);
n=n(NS_indep+1:end);
NS=NS-NS_indep;


% Compute the PDFs of the signals
[Px,X]=bs_sample_pdf_1d(x);
[Pn,N]=bs_sample_pdf_1d(n);
[Py,Y]=bs_sample_pdf_1d(y);

% Plot the PDF of each of the variables for inspection purpose.
clf;
subplot(3,1,1);
plot(X,Px,'linewidth',2,'color','b'); hold on
plot(Y,Py,'linewidth',2,'color','m');
plot(N,Pn,'linewidth',2,'color','g');
legend 'X' 'Y' 'N'

% Remove the sample mean from the random variables.
yp = y - sum(y)/NS;
xp = x - sum(x)/NS;

% Assume a linear relationship Y = alfa*X + N and
% compute alfa using the least square estimate LSQ.
alfa = sum(yp.*xp)/sum(xp.*xp);

% Now that you have alfa compute Yhat, you reconstruction
% of Y using the linear model.
yhat=xp*alfa;
% Add the sample mean of Y that was previously removed so that
% Yhat and Y can be compared.
yhat=yhat + sum(y)/NS;
% Compute the error associated with the estimate.
nhat= y - yhat; 
% Compute the sample PDF of the residuals and plot it.
[Pn,N]=bs_sample_pdf_1d( nhat );
plot(N,Pn,'linewidth',2,'color','k');
legend 'X' 'Y' 'N' 'Nhat'


% Compute the Joint PDF of X and Y.
[Pxy, X2, Y2]=bs_sample_pdf_2d(x,y);
% Plot the Joint PDF.
subplot(3,1,2); 
utl_contourfill(X2,Y2,Pxy); colorbar; hold on
title 'Joint PDF  Pxy(X,Y)'
xlabel 'X'
ylabel 'Y'

% Compare the LSQ estimate with the one inferred from the PDF.
%x1=input('Estimate Yhat. Select a value for X => ');
for iter=1:4
disp('Estimate Yhat. Select a value for X => ');
[x1,y1]=ginput(1);
% Perform estimate from JPDF.
y_jpdf=bs_estimate_pdf2_condx(Pxy, X2, Y2, x1);
y_jpdf
y1
% Perform LSQ estimate
y_lsq=alfa*x1;
% Plot the two estimates and check how close they are on the graph.
plot(x1,y_jpdf,'sb')
plot(x1,y_lsq,'*w');
disp('LSQ  estimate : black star')
disp('JPDF estimate : blue  square')
pause(0.2);
end

% Now use the JPDF and LSQ models to estimate Yhat with independent observations of Y
% that where not used to compute the JPDF and LSQ model.
subplot(3,1,3); 
clear y_lsq y_jpdf
for i=1:NS_indep;
	x1=x_indep(i);
    y_lsq(i) = alfa*x1;
    y_jpdf(i) = bs_estimate_pdf2_condx(Pxy, X2, Y2, x1);
end
% Make a plot 
hold on
plot(1:NS_indep, y_lsq,'linewidth',2,'color','r');
plot(1:NS_indep, y_jpdf,'linewidth',2,'color','g');
plot(1:NS_indep, y_indep,'color','b'); 

% Compute the skill associated with each estimate.
% Define skill = 1 - var(estimate - true)/ var(true)
mean_lsq=var(y_lsq(:)-y_indep(:)) ;
mean_jpdf=var(y_jpdf(:)-y_indep(:)) ;
skill_lsq  = 1 - mean_lsq /var( y_indep(:) );
skill_jpdf = 1 - mean_jpdf/var( y_indep(:) );
legend ( ['LSQ skill = ',num2str(skill_lsq)], ['Cond. JPDF = ',num2str(skill_jpdf)] ,'OBS');
xlabel 'X'
ylabel 'Y'
title 'Reconstruction of yhat'