% Optimal linear predictor for a correlated pair
% Textbook: Part I, "Linear Prediction", two-dimensional linear prediction example
%
% X = W1 and Y = 0.9*W1 + 0.1*W2, with W1, W2 ~ U[0,1] independent.
% The optimal linear predictor is Yhat = E[Y] + Cov(X,Y)/Var[X]*(x - E[X])
% = 1/2 + 0.9*(x - 1/2), and its MSE is Var[Y](1 - rho^2) = 0.01/12.

clearvars; close all
rng(0)

N = 1000;
W1 = rand(1, N);
W2 = rand(1, N);
X = W1;
Y = 0.9*W1 + 0.1*W2;

% Moments estimated from the data
E1 = mean(X);
E2 = mean(Y);
C12 = mean(X.*Y) - E1*E2;
V1 = var(X);
V2 = var(Y);
rho = C12/sqrt(V1*V2);

Yh = E2 + C12/V1*(X - E1);      % linear prediction
mse = mean((Y - Yh).^2);

figure
plot(X, Y, '.')
hold on
x = [0 1];
plot(x, E2 + C12/V1*(x - E1), 'r', 'LineWidth', 1.5)
hold off
axis equal
xlim([0 1])
xlabel('X')
ylabel('Y')
legend('data', 'linear predictor', 'Location', 'northwest')
grid on

fprintf('slope: theory = 0.9,    estimate = %.4f\n', C12/V1)
fprintf('rho:   theory = %.4f, estimate = %.4f\n', 0.9/sqrt(0.82), rho)
fprintf('MSE:   theory = %.3e, simulation = %.3e\n', 0.01/12, mse)