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

display(Copula_parameter);
%
% Number of copula parameters, Parameter_Copula
Parameter_Copula = 1;
cPlackett = @(u,v,theta)((1 + (theta-1)*(u + v))^2-4*u*v*theta*(theta-1))
^(-3/2)*theta*(1 + (theta-1)*(u + v-2*u*v));
c16 = @(u,v,theta)0.5*(1 + theta/u^2)*(1 + theta/v^2)*((u + v-1-theta*
(1/u + 1/v-1))^2 + 4*theta)^(-0.5)*(-((u + v-1-theta*(1/u + 1/v-1))^2 +
4*theta)^(-1)*(u + v-1-theta*(1/u + 1/v-1))^2 + 1);
%
% Calculate the AIC and BIC values for the four candidate copulas
PDF_Gaussian = copulapdf('Gaussian',U,rho_Gaussian(1,2));
AIC_Gaussian = -2*sum(log(PDF_Gaussian)) + 2*Parameter_Copula;
BIC_Gaussian = -2*sum(log(PDF_Gaussian)) + log(rows)*Parameter_Copula;
%
PDF_Plackett = zeros(rows,1);
for m = 1:rows
PDF_Plackett(m) = cPlackett(U(m,1),U(m,2),theta_Plackett);
end
AIC_Plackett = -2*sum(log(PDF_Plackett)) + 2*Parameter_Copula;
BIC_Plackett = -2*sum(log(PDF_Plackett)) + log(rows)*Parameter_Copula;
%
PDF_Frank = copulapdf('Frank',U,theta_Frank);
AIC_Frank = -2*sum(log(PDF_Frank)) + 2*Parameter_Copula;
BIC_Frank = -2*sum(log(PDF_Frank)) + log(rows)*Parameter_Copula;
%
PDF_16 = zeros(rows,1);
for m = 1:rows
PDF_16(m) = c16(U(m,1),U(m,2),theta_16);
end
AIC_16 = -2*sum(log(PDF_16)) + 2*Parameter_Copula;
BIC_16 = -2*sum(log(PDF_16)) + log(rows)*Parameter_Copula;
%
AIC_Copula = [AIC_Gaussian AIC_Plackett AIC_Frank AIC_16];
display(AIC_Copula);
[AIC_min_Copula Index_Copula] = min(AIC_Copula,[],2);
if Index_Copula = = 1
display('The best-fit copula using AIC is Gaussian copula');
elseif Index_Copula = = 2
display('The best-fit copula using AIC is Plackett copula');
elseif Index_Copula = = 3
display('The best-fit copula using AIC is Frank copula');
elseif Index_Copula = = 4
display('The best-fit copula using AIC is No.16 copula');
end
%
BIC_Copula = [BIC_Gaussian BIC_Plackett BIC_Frank BIC_16];
display(BIC_Copula);
[BIC_min_Copula Index_Copula] = min(BIC_Copula,[],2);
if Index_Copula = = 1
display('The best-fit copula using BIC is Gaussian copula');
elseif Index_Copula = = 2
display('The best-fit copula using BIC is Plackett copula');
elseif Index_Copula = = 3
display('The best-fit copula using BIC is Frank copula');