function plotSchematicFitting(ctrl,XInd,XP,x1Mean,x1Sigma,x2Mean,x2Sigma)

rng(ctrl.rngSeed);

if nargin < 6
    disp('Only formating inputs given. Using stored data.') 

    n = 5000;    % Number of samples to draw
    
    mu1 = 0;
    mu2 = 0;
    sigma1 = 0.7;        % covariance of margin 1
    sigma2 = 1.2;        % covariance of margin 2
    rho = -0.5;          % linear correlation coefficient \rho_{12}

    SigmaInd = [sigma1.^2 rho*sigma1*sigma2 ; rho*sigma1*sigma2 sigma2.^2];

    XInd = mvnrnd([mu1 mu2],SigmaInd,n);    % Draw from 
    XP = [1.2 -1];
end

[x1Mean,x1Sigma] = normfit(XInd(:,1));      % Prep, subplot (2,1)
[x2Mean,x2Sigma] = normfit(XInd(:,2));      % Prep, subplot (2,2)

% u1 = normcdf(XInd(:,1),x1Mean,x1Sigma);     % Prep, subplot (3,1)
% u2 = normcdf(XInd(:,2),x2Mean,x2Sigma);     % Prep, subplot (4,1)
% uXP1 = normcdf(XP(1),x1Mean,x1Sigma);       % Prep, subplot (3,1)
% uXP2 = normcdf(XP(2),x2Mean,x2Sigma);       % Prep, subplot (4,1)
u1 = normcdf(XInd(:,1),x1Mean,x1Sigma);     % Prep, subplot (3,1)
u2 = normcdf(XInd(:,2),x2Mean,x2Sigma);     % Prep, subplot (4,1)
uXP1 = normcdf(XP(1),x1Mean,x1Sigma);       % Prep, subplot (3,1)
uXP2 = normcdf(XP(2),x2Mean,x2Sigma);       % Prep, subplot (4,1)

xVec = -4:0.2:4;                            % x-axis for subplots in real space
u = linspace(0,1,100);                      % x-axis for subplot in uniform space
r1 = fitdist(XInd(:,1),'Normal');
r2 = fitdist(XInd(:,2),'Normal');

RhoEst = copulafit('Gaussian',[u1 u2],'Method','ML');   % Preparation of subplot (3:4,2)
[u1X,u2X] = meshgrid(u,u);
y = copulacdf('Gaussian',[u1X(:),u2X(:)],RhoEst);



figure('units','centimeters','OuterPosition',[10 10 20 20]);
tiledlayout(4, 2,'Padding','none');


% SUBPLOT (1:2,1)
nexttile([2,1]) 
plot(XInd(:,1),XInd(:,2),'.k')
hold on
plot(XP(1),XP(2),'dw','markerfacecolor',ctrl.c1,'markersize',8,'linewidth',1)
axis([-4 4 -4 4])
legend('$\mathbf{X} = \left[ \matrix{ X_1 \cr X_2 \cr}\right] \sim N\left(\left[ \matrix{ 0 \cr 0 \cr}\right] , \left[ \matrix{ \sigma_1^2 & \rho\sigma_1\sigma_2 \cr \rho\sigma_1\sigma_2 & \sigma_2^2 \cr}\right] \right)$','location','northeast','interpreter',ctrl.interpreter)
xlabel('$x_1$','interpreter',ctrl.interpreter)
ylabel('$x_2$','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


% SUBPLOT (1,2)
nexttile;
[a,~] = histcounts(XInd(:,1),xVec,'normalization','pdf');
f1t = find(xVec > XP(1),1,'first');
patch([xVec(f1t) xVec(f1t+1) xVec(f1t+1) xVec(f1t) xVec(f1t)],[0 0 a(f1t) a(f1t) 0],'w','FaceColor',ctrl.c1,'edgeColor','none','handleVisibility','off') 
hold on
histogram(XInd(:,1),xVec,'normalization','pdf','displaystyle','bar','edgecolor',[0 0 0],'facecolor',[1 1 1],'facealpha',0);
plot(xVec,pdf(r1,xVec),'-k','linewidth',1)
xlabel('$x_1$','interpreter',ctrl.interpreter); ylabel('$f_1(x_1)$','interpreter',ctrl.interpreter)
legend('Normalized histogram',['Fit $X_1 \sim \hat{N}($',sprintf('%2.1f,%2.1f',r1.mu,r1.sigma),'$)$'],'location','northeast','interpreter',ctrl.interpreter)
xlim([-4 4]); ylim([0 0.85])
set(gca,'ticklabelinterpreter',ctrl.interpreter)


% SUBPLOT (2,2)
nexttile;
[a,~] = histcounts(XInd(:,2),xVec,'normalization','pdf');
f1t = find(xVec > XP(2),1,'first');
patch([xVec(f1t) xVec(f1t+1) xVec(f1t+1) xVec(f1t) xVec(f1t)],[0 0 a(f1t) a(f1t) 0],'w','FaceColor',ctrl.c1,'edgeColor','none','handleVisibility','off') 
hold on 
histogram(XInd(:,2),xVec,'normalization','pdf','displaystyle','bar','edgecolor',[0 0 0],'facecolor',[1 1 1],'facealpha',0);
plot(xVec,pdf(r2,xVec),'-k','linewidth',1)
xlabel('$x_2$','interpreter',ctrl.interpreter); ylabel('$f_2(x_2)$','interpreter',ctrl.interpreter);
legend('Normalized histogram',['Fit $X_2 \sim \hat{N}($',sprintf('%2.1f,%2.1f',r2.mu,r2.sigma),'$)$'],'location','northeast','interpreter',ctrl.interpreter)
xlim([-4 4]); ylim([0 0.85])
set(gca,'ticklabelinterpreter',ctrl.interpreter)


% SUBPLOT (3,1)
nexttile;
[a,~] = histcounts(XInd(:,1),xVec,'normalization','cdf');
f1t = find(xVec > XP(1),1,'first');
patch([xVec(f1t) xVec(f1t+1) xVec(f1t+1) xVec(f1t) xVec(f1t)],[0 0 a(f1t) a(f1t) 0],'w','FaceColor',ctrl.c1,'edgeColor','none','handleVisibility','off') 
hold on
histogram(XInd(:,1),xVec,'normalization','cdf','displaystyle','bar','edgecolor',[0 0 0],'facecolor',[1 1 1],'facealpha',0);
hold on
plot(xVec,cdf(r1,xVec),'-k','linewidth',1)
xlabel('$x_1$','interpreter',ctrl.interpreter); ylabel('$F_1(x_1)$','interpreter',ctrl.interpreter)
xlim([-4 4]); ylim([0 1])
set(gca,'ticklabelinterpreter',ctrl.interpreter)


% SUBPLOT (4,1)
nexttile;
histogram2(u1,u2,linspace(0,1,10),linspace(0,1,10),'normalization','cdf','displaystyle','tile', ...
           'edgecolor','none','handlevisibility','off')
text(0.02,0.10,'2D cumulative histogram of $(F_1(x_1),F_2(x_2))$','interpreter',ctrl.interpreter,'fontsize',8,'backgroundcolor','w')
hold on
plot(uXP1,uXP2,'dw','markerfacecolor',ctrl.c1,'markersize',8,'linewidth',1,'handlevisibility','off')
xlabel('$F_1(x_1)$','interpreter',ctrl.interpreter)
ylabel('$F_2(x_2)$','interpreter',ctrl.interpreter)
colorbar('ticklabelinterpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)



nexttile;
[a,~] = histcounts(XInd(:,2),xVec,'normalization','cdf');
f1t = find(xVec > XP(2),1,'first');
patch([xVec(f1t) xVec(f1t+1) xVec(f1t+1) xVec(f1t) xVec(f1t)],[0 0 a(f1t) a(f1t) 0],'w','FaceColor',ctrl.c1,'edgeColor','none','handleVisibility','off') 
hold on
histogram(XInd(:,2),xVec,'normalization','cdf','displaystyle','bar','edgecolor',[0 0 0],'facecolor',[1 1 1],'facealpha',0);
plot(xVec,cdf(r2,xVec),'-k','linewidth',1)
xlabel('$x_2$','interpreter',ctrl.interpreter); ylabel('$F_2(x_2)$','interpreter',ctrl.interpreter)
xlim([-4 4]); ylim([0 1])
set(gca,'ticklabelinterpreter',ctrl.interpreter)



nexttile;
contour(u1X,u2X,reshape(y,100,100),'linewidth',1)
hold on
plot(uXP1,uXP2,'dw','markerfacecolor',ctrl.c1,'markersize',8,'linewidth',1,'handlevisibility','off')
view([0 90])
xlabel('$F_1(x_1)$','interpreter',ctrl.interpreter); ylabel('$F_2(x_2)$','interpreter',ctrl.interpreter)
legend(['Fit $C_{12}\left(F_1(x_1),F_2(x_2) ; \left[ \matrix{ ',sprintf('%2.1f & %2.1f',RhoEst(1),RhoEst(2)),' \cr ',sprintf('%2.1f & %2.1f',RhoEst(3),RhoEst(4)), ' \cr}\right] \right)$'], ...
       'location','southwest','interpreter',ctrl.interpreter)
colorbar('ticklabelinterpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


% Arrows between the different figures
annotation('arrow',[0.475 0.525],[0.75 0.75])
annotation('arrow',[0.525 0.475],[0.525 0.475])
annotation('arrow',[0.475 0.525],[0.30 0.30])


% Export
set(gcf,'PaperPositionMode','auto')
for aLoop = 1:numel(ctrl.plotFormat)
    print([ctrl.outputDirectory filesep 'figureFour'],ctrl.plotFormat{aLoop})
end