%% Script for blog post on correlation and partial correlation

%% 1. Generate some random variables of size N=250

nbins = 15;
fsize = 30;

f = 5.5;

m = 150;
s = 50;

N = 250;
X = m + randn(N,1) * s;
Y00 = m + randn(N,1) * s;
Yn10 = -(X - min(X)) + max(X); % invert
Yn05 = Yn10 + (f*(rand(N,1)-1)) * s;
Yn05 = (Yn05 - min(Yn05)) * (300 / range(Yn05));
Yp10 = X;
Yp05 = X + (f*(rand(N,1)-1)) * s;
Yp05 = (Yp05 - min(Yp05)) * (300 / range(Yp05));

%% 2. Plot with histograms

close all;

% r = -1
h = figure;
h.Color = 'w';
hh = scatterhist(X, Yn10, 'NBins', nbins);
hh = title('$r_{xy}=-1$', 'interpreter', 'latex');
hh.FontSize = fsize;
hh = xlabel('Variable X');
hh.FontSize = 16;
hh = ylabel('Variable Y');
hh.FontSize = 16;

% r = -0.5
h = figure;
h.Color = 'w';
hh = scatterhist(X, Yn05, 'NBins', nbins);
hh = title(sprintf('$r_{xy}=%1.1f$', corr(X,Yn05)), 'interpreter', 'latex');
hh.FontSize = fsize;
hh = xlabel('Variable X');
hh.FontSize = 16;
hh = ylabel('Variable Y');
hh.FontSize = 16;

% r = 0
h = figure;
h.Color = 'w';
hh = scatterhist(X, Y00, 'NBins', nbins);
hh = title('$r_{xy}=0.0$', 'interpreter', 'latex');
hh.FontSize = fsize;
hh = xlabel('Variable X');
hh.FontSize = 16;
hh = ylabel('Variable Y');
hh.FontSize = 16;

% r = 0.5
h = figure;
h.Color = 'w';
hh = scatterhist(X, Yp05, 'NBins', nbins);
hh = title(sprintf('$r_{xy}=%1.1f$', corr(X,Yp05)), 'interpreter', 'latex');
hh.FontSize = fsize;
hh = xlabel('Variable X');
hh.FontSize = 16;
hh = ylabel('Variable Y');
hh.FontSize = 16;

% r = 1
h = figure;
h.Color = 'w';
hh = scatterhist(X, X, 'NBins', nbins);
hh = title(sprintf('$r_{xy}=%1.1f$', corr(X,X)), 'interpreter', 'latex');
hh.FontSize = fsize;
hh = xlabel('Variable X');
hh.FontSize = 16;
hh = ylabel('Variable Y');
hh.FontSize = 16;


%% 3. Z-scores

close all;

clr_pos = [.2 .9 .2 .5];
clr_neg = [1 .9 .4 .5];

ntri = 3;
% idx = round(rand(ntri,1)*N);
idx = [46 109 225];

X_Z = zscore(X);
Yp10_Z = zscore(Yp10);
Yp05_Z = zscore(Yp05);
Yn10_Z = zscore(Yn10);

% r = 0.5
h = figure;
h.Position(3:4) = [560 500];
h.Color = 'w';
hh = scatterhist(X_Z, Yp05_Z, 'NBins', nbins);

hold on;
% Draw rectangles
for i = 1 : length(idx)
    x = X_Z(idx(i));
    y = Yp05_Z(idx(i));
    x1 = min(0, x);
    y1 = min(0, y);
    if x*y > 0
        clr = clr_pos;
    else
        clr = clr_neg;
    end
    hh = rectangle('Position', [x1, y1, abs(x), abs(y)], ...
                   'FaceColor', clr);
               
    plot(x, y, '.', 'MarkerSize', 50, 'Color', [0 .6 0]);
end

% Draw axes

hh = line([0 0], [-5 5]);
hh.Color = 'r';
hh.LineWidth = 2;
hh = line([-5 5], [0 0]);
hh.Color = 'r';
hh.LineWidth = 2;

hh = title(sprintf('$r_{xy}=%1.1f$', corr(X,Yp05)), 'interpreter', 'latex');
hh.FontSize = fsize;
hh = xlabel('Variable X (Z-score)');
hh.FontSize = 16;
hh = ylabel('Variable Y (Z-score)');
hh.FontSize = 16;
xlim([-3 3]);
ylim([-3 3]);

% r = 1
h = figure;
h.Position(3:4) = [560 500];
h.Color = 'w';
hh = scatterhist(X_Z, X_Z, 'NBins', nbins);

hold on;

% Draw rectangles
for i = 1 : length(idx)
    x = X_Z(idx(i));
    y = x;
    x1 = min(0, x);
    y1 = min(0, y);
    hh = rectangle('Position', [x1, y1, abs(x), abs(y)], ...
                   'FaceColor', clr_pos);
               
    plot(x, y, '.', 'MarkerSize', 50, 'Color', [0 .6 0]);
end

hh = line([0 0], [-5 5]);
hh.Color = 'r';
hh.LineWidth = 2;
hh = line([-5 5], [0 0]);
hh.Color = 'r';
hh.LineWidth = 2;

hh = title(sprintf('$r_{xy}=%1.1f$', corr(X,X)), 'interpreter', 'latex');
hh.FontSize = fsize;
hh = xlabel('Variable X (Z-score)');
hh.FontSize = 16;
hh = ylabel('Variable Y (Z-score)');
hh.FontSize = 16;
xlim([-3 3]);
ylim([-3 3]);


% r = -1
h = figure;
h.Position(3:4) = [560 500];
h.Color = 'w';
hh = scatterhist(X_Z, Yn10_Z, 'NBins', nbins);

hold on;

% Draw rectangles
for i = 1 : length(idx)
    x = X_Z(idx(i));
    y = Yn10_Z(idx(i));
    x1 = min(0, x);
    y1 = min(0, y);
    hh = rectangle('Position', [x1, y1, abs(x), abs(y)], ...
                   'FaceColor', clr_neg);
               
    plot(x, y, '.', 'MarkerSize', 50, 'Color', [0 .6 0]);
end

hh = line([0 0], [-5 5]);
hh.Color = 'r';
hh.LineWidth = 2;
hh = line([-5 5], [0 0]);
hh.Color = 'r';
hh.LineWidth = 2;

hh = title(sprintf('$r_{xy}=%1.1f$', corr(X,Yn10)), 'interpreter', 'latex');
hh.FontSize = fsize;
hh = xlabel('Variable X (Z-score)');
hh.FontSize = 16;
hh = ylabel('Variable Y (Z-score)');
hh.FontSize = 16;
xlim([-3 3]);
ylim([-3 3]);

%% 4. Invert a correlation matrix to compute partial correlations

M=[[1 .43 .5];[.43 1 .39];[.5 .39 1]];
P=inv(M);

Mp = zeros(3);

for i = 1 : 3
   for j = 1 : 3
       if i ~= j
           Mp(i,j) = -P(i,j) / sqrt(P(i,i)*P(j,j));
           Mp(j,i) = Mp(i,j);
       end
   end
end

Mp
