Modeling and simulation of bivariate distribution of shear strength parameters using copulas - Risk and Reliability in Geotechnical Engineering - page 115

identification of the best-fit marginal distributions and copula underlying measured data,
(2) MATLAB codes for simulation of copulas and bivariate distribution, and (3) MATLAB
codes for definition of subfunctions used in (1) and (2).
(1) identification of best-fit margins and copula using aiC and BiC
% Identification of best-fit margins and copula using AIC and BIC
% Filename: Identification.m
% © (2013) Dian-Qing Li and Xiao-Song Tang
clc;
format short;
% Read observed cohesions and friction angles stored in XLS-files, data
data = xlsread('data.xls');
%
% Select best-fit margin among TruncNormal,Lognormal,TruncGumbel and
% Weibull distributions
% Rows,columns,mean,standard deviation of data, rows,cols,mu,sigma
[rows cols] = size(data);
mu = mean(data);
sigma = std(data);
% Number of distribution parameters, Parameter_Margin
Parameter_Margin = 2;
%
% Calculate the AIC and BIC values for variable X1
PDF_TruncNormal_X1 = normpdf(data(:,1),mu(1),sigma(1))/(1-normcdf(-mu(1)/
sigma(1)));
AIC_TruncNormal_X1 = -2*sum(log(PDF_TruncNormal_X1)) + 2*Parameter_Margin;
BIC_TruncNormal_X1 = -2*sum(log(PDF_TruncNormal_X1)) + log(rows)*Parameter_
Margin;
%
sLn = sqrt(log(1 + (sigma(1)/mu(1))^2)); mLn = log(mu(1))-sLn^2/2;
PDF_Lognormal_X1 = lognpdf(data(:,1),mLn,sLn);
AIC_Lognormal_X1 = -2*sum(log(PDF_Lognormal_X1)) + 2*Parameter_Margin;
BIC_Lognormal_X1 = -2*sum(log(PDF_Lognormal_X1)) + log(rows)*Parameter_
Margin;
%
a = 1.282549808/sigma(1); b = mu(1)-(0.5772156649/a);
PDF_TruncGumbel_X1 = a*exp(-a*(data(:,1)-b)-exp(-a*(data(:,1)-b)))/
(1-exp(-exp(a*b)));
AIC_TruncGumbel_X1 = -2*sum(log(PDF_TruncGumbel_X1)) + 2*Parameter_Margin;
BIC_TruncGumbel_X1 = -2*sum(log(PDF_TruncGumbel_X1)) + log(rows)*Parameter_
Margin;
%
k = k_solveWeibull(mu(1),sigma(1)); u = mu(1)/gamma(1 + 1/k);
PDF_Weibull_X1 = k/u*(data(:,1)/u).^(k-1).*exp(-(data(:,1)/u).^k);
AIC_Weibull_X1 = -2*sum(log(PDF_Weibull_X1)) + 2*Parameter_Margin;
BIC_Weibull_X1 = -2*sum(log(PDF_Weibull_X1)) + log(rows)*Parameter_Margin;
%
AIC_X1 = [AIC_TruncNormal_X1 AIC_Lognormal_X1 AIC_TruncGumbel_X1 AIC_
Weibull_X1];
display(AIC_X1);
[AIC_min_X1 Index_X1] = min(AIC_X1,[],2);
if Index_X1 = = 1
display('The best-fit margin for X1 using AIC is TruncNormal
distribution');