
% ----------------------------------------------------
% explanation
% answers to Q7 for Gibbs sampling
% ----------------------------------------------------

%% housekeeping

clc;clear;close all;
set(0,'defaulttextinterpreter','latex');


%% Parameter

% mean 
mu = [1;2];
% variance-covariance matrix
Sigma = rand(2)/10;
% make it positive definite
Sigma = Sigma'*Sigma;

% number of simulations
num_sim = 20000;

% box. This is called M in the pdf.
X = nan(2,num_sim);

% initial values (irrelevant for large num_sim)
X(:,1) = mu;

%% Gibbs sampling part

for i = 2 : num_sim
	
    % conditional distribution for variable 1 given variable 2
    m_x1_given_x2 = mu(1) + Sigma(1,2)/Sigma(2,2) * (X(2,i-1)-mu(2));
    sigma_x1_given_x2 = (Sigma(1,1)*Sigma(2,2)-Sigma(1,2)^2)/Sigma(2,2);
    X(1,i) =m_x1_given_x2 + sqrt(sigma_x1_given_x2) * randn;
    
    % conditional distribution for variable 2 given variable 1
    m_x2_given_x1 = mu(2) + Sigma(1,2)/Sigma(1,1) * (X(1,i)-mu(1));
    sigma_x2_given_x1 = (Sigma(1,1)*Sigma(2,2)-Sigma(1,2)^2)/Sigma(1,1);
    X(2,i) = m_x2_given_x1 + sqrt(sigma_x2_given_x1) * randn;
    
    % small checks.
    if abs(X(1,i)) > 1e10 || abs(X(2,i)) > 1e10
    	keyboard;
    end

end

%% True densities

x1 = min(X(1,:)) : .001: max(X(1,:));
phi_x1 = normpdf(x1,mu(1),sqrt(Sigma(1,1)));

x2 = min(X(2,:)) : .001: max(X(2,:));
phi_x2 = normpdf(x2,mu(2),sqrt(Sigma(2,2)));

%% Graph
figure;

subplot(2,1,1)
hold on
h1 = histogram(X(1,:));
h1.Normalization = 'pdf';
plot(x1,phi_x1,'-r','LineWidth',2.5);
legend('Histogram','True density','Location','Best')
title('$X_1$')

subplot(2,1,2)
hold on
h1 = histogram(X(2,:));
h1.Normalization = 'pdf';
plot(x2,phi_x2,'-r','LineWidth',2.5);
legend('Histogram','True density','Location','Best')
title('$X_2$')

%% Correlation
disp('Simulated Cov(X1,X2)=')
disp(cov(X'))
disp('True Cov(X1,X2) ( = Sigma ) =')
disp(Sigma)
1+1;
