% Gaussian marginals do not make a jointly Gaussian pair
% Textbook: Part I, "Two Random Variables", the bivariate Gaussian
%
% X ~ N(0,1) and W = +-1 equally likely, independent of X. Then Y = W*X is
% also N(0,1), and Cov(X,Y) = 0. Yet X and Y are dependent (|Y| = |X|), and
% the pair is not jointly Gaussian: Z = X + Y equals 0 with probability 1/2,
% which no Gaussian variable can do.

clearvars; close all
rng(0)

N = 1e6;
X = randn(1, N);
W = 2*randi(2, [1 N]) - 3;      % +-1 values
Y = W.*X;
Z = X + Y;

figure
tiledlayout(2, 2)
nexttile; histogram(X, 'Normalization', 'pdf'); title('X'); grid on
nexttile; histogram(Y, 'Normalization', 'pdf'); title('Y = WX'); grid on
nexttile; plot(X(1:2000), Y(1:2000), '.'); axis equal; xlabel('X'); ylabel('Y'); grid on
title('joint sample')
nexttile; histogram(Z, 'Normalization', 'probability'); title('Z = X + Y'); grid on

fprintf('Var[Y]:   theory = 1, simulation = %.4f\n', var(Y))
fprintf('Cov(X,Y): theory = 0, simulation = %.4f\n', mean(X.*Y) - mean(X)*mean(Y))
fprintf('Pr(Z = 0): %.4f (not Gaussian)\n', mean(Z == 0))