% Linear prediction for a bivariate Gaussian pair
% Textbook: Part I, "Linear Prediction"
%
% (X,Y) jointly Gaussian with covariance matrix C. The predictor
% Yhat = E[Y] + Cov(X,Y)/Var[X]*(x - E[X]) has MSE = Var[Y](1 - rho^2).
% For a Gaussian pair this linear predictor is also the best predictor of
% any form.

clearvars; close all
rng(0)

mu = [0; 0];
C = [2 0.75; 0.75 1];       % Var[X] = 2, Var[Y] = 1, Cov(X,Y) = 0.75
N = 1e4;

rho = C(1,2)/sqrt(C(1,1)*C(2,2));

L = chol(C, 'lower');       % correlated pairs (or mvnrnd(mu', C, N))
Z = mu + L*randn(2, N);
X = Z(1,:);
Y = Z(2,:);

Yh = mu(2) + C(1,2)/C(1,1)*(X - mu(1));
e = Y - Yh;

figure
plot(X, Y, '.')
hold on
x = [-5 5];
plot(x, mu(2) + C(1,2)/C(1,1)*(x - mu(1)), 'r', 'LineWidth', 1.5)
hold off
axis equal
xlabel('X')
ylabel('Y')
legend('data', 'linear predictor', 'Location', 'northwest')
grid on

fprintf('rho = %.4f\n', rho)
fprintf('MSE: theory = %.4f, simulation = %.4f\n', (1 - rho^2)*C(2,2), mean(e.^2))