% Joint CDF of a bivariate Gaussian
% Textbook: Part I, "Two Random Variables", the bivariate Gaussian
%
% F(a,b) = Pr(X <= a, Y <= b) for (X,Y) jointly Gaussian with means mu,
% standard deviations sigma1, sigma2 and correlation coefficient rho.
% Estimated by simulation, and compared with mvncdf (Statistics and Machine
% Learning Toolbox). At the mean, F(mu1,mu2) = 1/4 + asin(rho)/(2*pi).

clearvars; close all
rng(0)

mu = [0; 0];
sigma1 = 1;
sigma2 = 2;
rho = 0.5;
C = [sigma1^2, rho*sigma1*sigma2; rho*sigma1*sigma2, sigma2^2];  % covariance matrix

% Correlated Gaussian pairs from independent N(0,1) (or mvnrnd(mu', C, N))
N = 1e6;
L = chol(C, 'lower');
Z = mu + L*randn(2, N);

a = 0; b = 0;
F_sim = mean(Z(1,:) <= a & Z(2,:) <= b);
fprintf('F(%g,%g): simulation = %.4f\n', a, b, F_sim)
fprintf('F(%g,%g): 1/4 + asin(rho)/(2 pi) = %.4f\n', a, b, 1/4 + asin(rho)/(2*pi))
if exist('mvncdf', 'file')
    fprintf('F(%g,%g): mvncdf = %.4f\n', a, b, mvncdf([a b], mu', C))
end

figure
histogram2(Z(1,:), Z(2,:), 'Normalization', 'pdf', 'DisplayStyle', 'tile')
axis equal
xlabel('x')
ylabel('y')
title('f_{X,Y}(x,y), \rho = 0.5')
colorbar