MATLAB中Copula函数拟合代码源码分享

代码简介：

数据并非来自个人，故不公开，0基础需跑码者，请自学copula函数相关理论、概率论中关于边缘分布和联合概率分布的相关知识、MATLAB中内置的边缘概率分布拟合等操作基础。本代码中设置变量过多，可能略有阅读难度，敬请谅解。

本伪代码默认已得到边缘概率分布拟合结果：本伪代码使用的变量j-Precipitation降水变量分布为Lognormal、变量q-Temperature气温分布为burr，在此基础上取样100个点生成二者的cdf值，作为copula函数的自变量，对GussianCopula、tCopula、GumbelCopula、ClaytonCopula、FrankCopula共5个MATLAB内置函数使用copulafit功能做出拟合，并绘制Copula函数密度图、分布图以及自变量为降水和气温的二元联合概率密度图、分布图，并没有对拟合的copula函数做出拟合优度检验（请参考谢中华老师编写的《MATLAB统计分析与应用：40个案例分析》，本人也主要参考其中代码学习进行copula拟合）

代码部分：

分为两个.m文件，第一个文件为数据处理生成必要变量的文件，以第二个copula文件为主，且其中部分变量在第一个文件中声明和赋值

①数据处理及基础变量文件

UY0225 = readmatrix('UY0225.xlsx');

j = UY0225(:,2);

q = UY0225(:,3);

j_ksdensity_pdf = ksdensity(j,j,"Function","pdf");

q_ksdensity_pdf = ksdensity(q,q,"Function","pdf");

j_ksdensity_cdf = ksdensity(j,j,'Function','cdf');

q_ksdensity_cdf = ksdensity(q,q,'Function','cdf');

[j_ecdf,jsort_ecdf] = ecdf(j);%计算经验分布函数值

[q_ecdf,qsort_ecdf] = ecdf(q);

j_ecdf_spline = spline(jsort_ecdf(2:end),j_ecdf(2:end),j);%对ecdf做样本点差值

q_ecdf_spline = spline(qsort_ecdf(2:end),q_ecdf(2:end),q);

j_axis = linspace(min(j),max(j),100);%%%100—50

q_axis = linspace(min(q),max(q),100);%%%100-50

[j_grid,q_grid] = meshgrid(j_axis,q_axis);

j_ecdf_spline_axis = spline(jsort_ecdf(2:end),j_ecdf(2:end),j_axis);%对ecdf做坐标轴差值

q_ecdf_spline_axis = spline(qsort_ecdf(2:end),q_ecdf(2:end),q_axis);

[j_ecdf_spline_grid,q_ecdf_spline_grid] = meshgrid(j_ecdf_spline_axis,q_ecdf_spline_axis);

[jsort,idj] = sort(j);%对j排序

[qsort,idq] = sort(q);

j_axis_ksdensity = ksdensity(j,j_axis,"Function","cdf");%核密度计算j_axis处的值

q_axis_ksdensity = ksdensity(q,q_axis,"Function","cdf");%pts位子内必须是向量点，不能是矩阵

[j_grid_ksdensity,q_grid_ksdensity] = meshgrid(j_axis_ksdensity,q_axis_ksdensity);

②copula函数拟合文件

eval(['basic_data_calculation;']);

%----------最佳mpdf----------

Log_normal_j = fitdist(j,'Lognormal');

bur_q = fitdist(q,'Burr');

j_Lognormal_axis = cdf(Log_normal_j,j_axis);

q_burr_axis = cdf(bur_q,q_axis);

[j_Lognormal_grid,q_burr_grid] = meshgrid(j_Lognormal_axis,q_burr_axis);

%二元联合概率分布做准备

j_bestmpdf_min = icdf(Log_normal_j,0.0001);%计算mpdf头尾值

j_bestmpdf_max = icdf(Log_normal_j,0.9999);

q_bestmpdf_min = icdf(bur_q,0.0001);

q_bestmpdf_max = icdf(bur_q,0.9999);

j_bestmpdf_jaxis = linspace(j_bestmpdf_min,j_bestmpdf_max,100);%生成坐标轴上等间格

q_bestmpdf_qaxis = linspace(q_bestmpdf_min,q_bestmpdf_max,100);

[j_bestmpdf_jgrid,q_bestmpdf_qgrid] = meshgrid(j_bestmpdf_jaxis,q_bestmpdf_qaxis);%生成网格

j_bestmpdf_jaxis_cdf = cdf(Log_normal_j,j_bestmpdf_jaxis);

q_bestmpdf_qaxis_cdf = cdf(bur_q,q_bestmpdf_qaxis);

[j_bestmpdf_jaxis_cdfgrid,q_bestmpdf_qaxis_cdfgrid] = meshgrid(j_bestmpdf_jaxis_cdf,q_bestmpdf_qaxis_cdf);

%----------拟合Guassian、t、Gumbel、Clayton、Frank的Copula's paramaters----------

rho_norm_LognormalBurr = copulafit('Gaussian',[j_Lognormal_axis(:),q_burr_axis(:)]);

[rho_t_LognormalBurr,nuci_LognormalBurr] = copulafit('t',[j_Lognormal_axis(:),q_burr_axis(:)]);

alpha_Gumbel_LognormalBurr = copulafit('Gumbel',[j_Lognormal_axis(:),q_burr_axis(:)]);

alpha_Clayton_LognormalBurr = copulafit('Clayton',[j_Lognormal_axis(:),q_burr_axis(:)]);

alpha_Frank_LognormalBurr = copulafit('Frank',[j_Lognormal_axis(:),q_burr_axis(:)]);

%----------计算copula(u,v)的值（以生成的等间距uv为自变量计算）----------

[j_LognormalBurr_grid,q_LognormalBurr_grid] = meshgrid(linspace(0,1));

Cpdf_norm_LognormalBurr = copulapdf("Gaussian",[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],rho_norm_LognormalBurr);

Ccdf_norm_LognormalBurr = copulacdf("Gaussian",[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],rho_norm_LognormalBurr);

Cpdf_t_LognormalBurr = copulapdf("t",[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],rho_t_LognormalBurr,nuci_LognormalBurr);

Ccdf_t_LognormalBurr = copulacdf("t",[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],rho_t_LognormalBurr,nuci_LognormalBurr);

Cpdf_Gumbel_LognormalBurr = copulapdf('Gumbel',[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],alpha_Gumbel_LognormalBurr);

Ccdf_Gumbel_LognormalBurr = copulacdf('Gumbel',[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],alpha_Gumbel_LognormalBurr);

Cpdf_Clayton_LognormalBurr = copulapdf('Clayton',[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],alpha_Clayton_LognormalBurr);

Ccdf_Clayton_LognormalBurr = copulacdf('Clayton',[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],alpha_Clayton_LognormalBurr);

Cpdf_Frank_LognormalBurr = copulapdf('Frank',[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],alpha_Frank_LognormalBurr);

Ccdf_Frank_LognormalBurr = copulacdf('Frank',[j_LognormalBurr_grid(:),q_LognormalBurr_grid(:)],alpha_Frank_LognormalBurr);

%----------计算F(j,p)的值（以生成的等间距j_bestmpdf,q_bestmpdf为自变量计算）----------

Cpdf_norm_LognormalBurr_grid = copulapdf("Gaussian",[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],rho_norm_LognormalBurr);

Ccdf_norm_LognormalBurr_grid = copulacdf('Gaussian',[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],rho_norm_LognormalBurr);

Cpdf_t_LognormalBurr_grid = copulapdf("t",[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],rho_t_LognormalBurr,nuci_LognormalBurr);

Ccdf_t_LognormalBurr_grid = copulacdf('t',[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],rho_t_LognormalBurr,nuci_LognormalBurr);

Cpdf_Gumbel_LognormalBurr_grid = copulapdf('Gumbel',[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],alpha_Gumbel_LognormalBurr);

Ccdf_Gumbel_LognormalBurr_grid = copulacdf('Gumbel',[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],alpha_Gumbel_LognormalBurr);

Cpdf_Clayton_LognormalBurr_grid = copulapdf('Clayton',[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],alpha_Clayton_LognormalBurr);

Ccdf_Clayton_LognormalBurr_grid = copulacdf('Clayton',[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],alpha_Clayton_LognormalBurr);

Cpdf_Frank_LognormalBurr_grid = copulapdf('Frank',[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],alpha_Frank_LognormalBurr);

Ccdf_Frank_LognormalBurr_grid = copulacdf('Frank',[j_bestmpdf_jaxis_cdfgrid(:),q_bestmpdf_qaxis_cdfgrid(:)],alpha_Frank_LognormalBurr);

%绘图

%--------正态copula和正态二元联合--------

figure;

subplot(2,2,1);%正态密度

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Cpdf_norm_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('c_n_o_r_m(u,v)');

title("c(u,v)——normCopula");

subplot(2,2,2);%正态分布

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Ccdf_norm_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('C_n_o_r_m(u,v)');

title("C(u,v)——normCopula");

subplot(2,2,3);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Cpdf_norm_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('f_n_o_r_m_a_l(P,T)');

title("f(P,T)——二元normCopula");

subplot(2,2,4);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Ccdf_norm_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('F_n_o_r_m_a_l(P,T)');

title("F(P,T)——二元normCopula");

%--------t-copula和t联合概率分布

figure;

subplot(2,2,1);%t密度

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Cpdf_t_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('c_t(u,v)');

subplot(2,2,2);%t分布

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Ccdf_t_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('C_t(u,v)');

subplot(2,2,3);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Cpdf_t_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('f_t(P,T)');

title("f(P,T)——二元tCopula");

subplot(2,2,4);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Ccdf_t_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('F_t(P,T)');

title("F(P,T)——二元tCopula");

%--------GumbelCopula和Gumbel二元联合--------

figure;

subplot(2,2,1);%GumbelCopula密度

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Cpdf_Gumbel_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('c_G_u_m_b_e_l(u,v)');

title("c(u,v)——GumbelCopula");

subplot(2,2,2);%GumbelCopula分布

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Ccdf_Gumbel_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('C_G_u_m_b_e_l(u,v)');

title("C(u,v)——GumbelCopula");

subplot(2,2,3);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Cpdf_Gumbel_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('f_G_u_m_b_e_l_a_l(P,T)');

title("f(P,T)——二元GumbelCopula");

subplot(2,2,4);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Ccdf_Gumbel_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('F_G_u_m_b_e_l_a_l(P,T)');

title("F(P,T)——二元GumbelCopula");

%--------ClaytonCopula和Clayton二元联合--------

figure;

subplot(2,2,1);%ClaytonCopula密度

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Cpdf_Clayton_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('c_C_l_a_y_t_o_n(u,v)');

title("c(u,v)——ClaytonCopula");

subplot(2,2,2);%ClaytonCopula分布

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Ccdf_Clayton_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('C_C_l_a_y_t_o_n(u,v)');

title("C(u,v)——ClaytonCopula");

subplot(2,2,3);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Cpdf_Clayton_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('f_C_l_a_y_t_o_n_a_l(P,T)');

title("f(P,T)——二元ClaytonCopula");

subplot(2,2,4);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Ccdf_Clayton_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('F_C_l_a_y_t_o_n_a_l(P,T)');

title("F(P,T)——二元ClaytonCopula");

%--------FrankCopula和Frank二元联合--------

figure;

subplot(2,2,1);%FrankCopula密度

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Cpdf_Frank_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('c_F_r_a_n_k(u,v)');

title("c(u,v)——FrankCopula");

subplot(2,2,2);%FrankCopula分布

surf(j_LognormalBurr_grid,q_LognormalBurr_grid,reshape(Ccdf_Frank_LognormalBurr,size(j_LognormalBurr_grid)));

xlabel('U(precipitation)');

ylabel('V(temperature)');

zlabel('C_F_r_a_n_k(u,v)');

title("C(u,v)——FrankCopula");

subplot(2,2,3);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Cpdf_Frank_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('f_F_r_a_n_k_a_l(P,T)');

title("f(P,T)——二元FrankCopula");

subplot(2,2,4);

surf(j_bestmpdf_jgrid,q_bestmpdf_qgrid,reshape(Ccdf_Frank_LognormalBurr_grid,size(j_bestmpdf_jgrid)));

xlabel('Precipitation');

ylabel('Temperature');

zlabel('F_F_r_a_n_k_a_l(P,T)');

title("F(P,T)——二元FrankCopula");

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删