%% Infer Marginals and Copulas

%% 1 - INITIALIZE UQLAB
clearvars
rng(100,'twister')
uqlab

%% 2 - DATA GENERATION
% Specify the marginals of these random variables:
iOptsTrue.Marginals(1).Name = 'X1';
iOptsTrue.Marginals(1).Type = 'Uniform';
iOptsTrue.Marginals(1).Moments = [0, 28.0];

iOptsTrue.Marginals(2).Name = 'X2';
iOptsTrue.Marginals(2).Type = 'Uniform';
iOptsTrue.Marginals(2).Moments = [0, 1];

iOptsTrue.Marginals(3).Name = 'X3';
iOptsTrue.Marginals(3).Type = 'Lognormal';
iOptsTrue.Marginals(3).Moments = [82.8 , 22.1];
iOptsTrue.Marginals(3).Bounds = [40.0 , 120];

iOptsTrue.Marginals(4).Name = 'X4';
iOptsTrue.Marginals(4).Type = 'Gaussian';
iOptsTrue.Marginals(4).Moments = [10.0, 6.0];
iOptsTrue.Marginals(4).Bounds = [1.0, 30.0];

iOptsTrue.Marginals(5).Name = 'X5';
iOptsTrue.Marginals(5).Type = 'Gaussian';
iOptsTrue.Marginals(5).Moments = [1.0, 0.5];
iOptsTrue.Marginals(5).Bounds = [0.1 , 2.0];

iOptsTrue.Marginals(6).Name = 'X6';
iOptsTrue.Marginals(6).Type = 'Gaussian';
iOptsTrue.Marginals(6).Moments = [2400, 100];
iOptsTrue.Marginals(6).Bounds = [2200 , 2600];

iOptsTrue.Marginals(7).Name = 'X7';
iOptsTrue.Marginals(7).Type = 'Gaussian';
iOptsTrue.Marginals(7).Moments = [1850, 50];
iOptsTrue.Marginals(7).Bounds = [1750 , 2000];

iOptsTrue.Marginals(8).Name = 'X8';
iOptsTrue.Marginals(8).Type = 'Lognormal';
iOptsTrue.Marginals(8).Moments = [0.62, 0.17];
iOptsTrue.Marginals(8).Bounds = [0.2, 2.0];

iOptsTrue.Marginals(9).Name = 'X9';
iOptsTrue.Marginals(9).Type = 'Gaussian';
iOptsTrue.Marginals(9).Moments = [30, 7];
iOptsTrue.Marginals(9).Bounds = [15 , 45];

iOptsTrue.Marginals(10).Name = 'X10';
iOptsTrue.Marginals(10).Type = 'Lognormal';
iOptsTrue.Marginals(10).Moments = [0.71, 0.13];
iOptsTrue.Marginals(10).Bounds = [0.5 , 2.0];

iOptsTrue.Marginals(11).Name = 'X11';
iOptsTrue.Marginals(11).Type = 'Gaussian';
iOptsTrue.Marginals(11).Moments = [34, 6];
iOptsTrue.Marginals(11).Bounds = [20.0, 45.0];

iOptsTrue.Marginals(12).Name = 'X12';
iOptsTrue.Marginals(12).Type = 'Uniform';
iOptsTrue.Marginals(12).Moments = [0.01, 0.99];

iOptsTrue.Marginals(13).Name = 'X13';
iOptsTrue.Marginals(13).Type = 'Gaussian';
iOptsTrue.Marginals(13).Moments = [30, 8];
iOptsTrue.Marginals(13).Bounds = [15.0 , 45.0];

iOptsTrue.Marginals(14).Name = 'X14';
iOptsTrue.Marginals(14).Type = 'Uniform';
iOptsTrue.Marginals(14).Moments = [100, 150];

iOptsTrue.Marginals(15).Name = 'X15';
iOptsTrue.Marginals(15).Type = 'Gaussian';
iOptsTrue.Marginals(15).Moments = [300, 150];
iOptsTrue.Marginals(15).Bounds = [100 ,1500];

iOptsTrue.Marginals(16).Name = 'X16';
iOptsTrue.Marginals(16).Type = 'Gaussian';
iOptsTrue.Marginals(16).Moments = [0.4, 0.3];
iOptsTrue.Marginals(16).Bounds = [0.01 1.0];

%--------------------------Independent Copula--------------------------
iOptsTrue.Copula.Type = 'Independent';

%%
% Create an INPUT object based on the specified marginals and copulas:
iOptsTrue.Name = 'True Input';
myInputTrue = uq_createInput(iOptsTrue);

%%
% Visualize the input model:
uq_print(myInputTrue)
%uq_figure
uq_display(myInputTrue)

%%
% Sample from the input model: 
X = uq_getSample(myInputTrue,400);

%%
U = uq_all_cdf(X,myInputTrue.Marginals);

%% 3 ---------------------------------- ok KS ok BIC
% clear iOpts
%-----------------------------------------------True Input is finished
iOpts.Inference.Data = X;
%--------------------------Kolmogorov-Smirnov (KS) -Marginals--------------------------
iOpts.Inference.Criterion = 'KS';
%--------------------------Bayesian (BIS) selection -Copula--------------------------
iOpts.Copula.Inference.Criterion = 'BIC';
iOpts.Name = 'CopulaForm';
CopulaForm = uq_createInput(iOpts);

%%
% Print out a report on the inferred input model:
uq_print(CopulaForm)
uq_display(CopulaForm)


%% 4 - VALIDATION OF BLOCK INDEPENDENCE
Blocks = {CopulaForm.Copula.Variables};
NrBlocks = length(Blocks);
VarNames = {CopulaForm.Marginals.Name};

%%
% Draw a sample from the input model:
XX = uq_getSample(1000);

%%
% Calculate the correlation matrix |R| on the sample |X| obtained above
R = corr(XX, 'type', 'Kendall');

fprintf('Full correlation matrix:\n')
uq_printMatrix(R,VarNames,VarNames)

%% 5 - STATISTICAL TEST OF BLOCK INDEPENDENCE
BlocksHat = uq_test_block_independence(XX, 0.05);



