(1) 多变量联合分布函数Copula 二维copula重现期matlab制图代码

1 多变量联合分布函数Copula

Copula函数是一个多变量统计方法，它允许研究和描述多个随机变量之间的依赖关系。该方法适用于什么情景？当有多个变量且变量不独立时，往往需要多个变量联合去表征某种信息，即为复合事件，Copula在水文水资源领域应用广泛，诸如洪水事件属性包括洪峰、洪量和洪水历时等，干旱事件属性包括干旱历时、干旱面积和干旱强度等，Copula函数可用于计算这些干旱特征或者洪水特征的联合和条件概率分布和重现期。

2 原理解析及计算流程

Copula理论可以追溯到1959年，Sklar通过定理形式提出，将一个联合分布分解为n个边缘分布和一个Copula函数，此时Copula函数可以描述这n个变量间的相关结构，构造灵活且求解相对简单。Copula函数是求解多变量概率问题优良的数学工具，通过独立随机变量的边缘分布来反映变量之间的相关性。Copula函数完全适用于多因素耦合安全风险分析，构建的多变量概率模型简便、适用性强，可以准确、客观地评估耦合安全因子协同作用下的发生风险概率大小，是减缓或扭转风险持续增长趋势的前提基础，近年来常被应用于多维致灾因子风险评估当中，如对洪水、干旱等事件进行风险研究。

整体来看，开展Copula函数应用的基本流程包括：①根据实际研究选定随机变量；②确定各随机变量的边缘分布；③计算相关系数并进行显著性检验；④确定备选Copula函数的参数；⑤对Copula函数进行拟合检验与优选；⑥推求联合概率分布和条件概率分布。

以干旱事件为例，在本研究中关键步骤主要为：

（1）构建干旱历时(H)、干旱面积(A)和干旱强度(S)的单变量边缘分布函数，开展显著性检验后确定最优拟合函数并计算累积概率；

（2）构建H-A、H-S、A-S的二维及H-A-S的三维联合分布函数，开展显著性检验后确定最优拟合函数并计算概率；

（3）通过RMSE、NSE、AIC、BIC等指标检验拟合效果，并得到最终概率计算。

计算流程图如图1所示。图1 流程图

3 技术名词解释

（1）单变量边缘分布

Copula函数不限定变量的边缘分布，通过Copula模型，可以将k个任意的边际分布连接起来，形成一个多变量联合分布概率模型。Copula函数的理论基础是Sklar’s定理，随机变量X1、X2、X3、…Xn连续，F1（x1）、F2（x2）、F3（x3）、…、Fn（xn）是其边缘分布函数，F为n维联合概率分布函数，则存在Copula函数C:[0,1]n→[0,1]，使得：

若边缘分布函数F1（x1）、F2（x2）、F3（x3）、…、Fn（xn）是连续的，则Copula函数唯一确定。

本文采用常见的单变量分布函数，如Weibull-Copula、Gamma-Copula、Exp-Copula、Normal-Copula和Lognormal-Copula等，采用极大似然法估算参数，并通过拟合优度检验（goodness-o-fit test）来判断一组数据是否服从于以上某种分布，常用的拟合优度检验方法有Q-Q图法、Kolmogorov-Smimov(K-S)检验等，其中K-S检验(Kolmogorov-Smirnov test)结果较为精准，本文采用K-S检验确立最优的边缘分布函数。各分布函数计算公式如表1所示：

表1 单变量边缘分布常见的Copula函数

K-S检验统计量D的定义如下：

式中ck为联合观测值样本（xk，yk）的Copula值；mk为联合观测值样本中满足条件x≤xk且y≤yk的联合观测值的个数。

（2）多变量联合分布

Copula作为联合分布函数和边缘分布函数之间的连接函数，包含了很多分布族，其中椭圆分布族和Archimedean分布族是常见的两个分布族，Archimedean分布族在水资源领域应用较多，其主要包括Frank-Copula、Clayton-Copula和Gumbel-Copula，二维及三维计算公式如表2所示。
表2 多变量联合分布常见的Archimedean Copula函数

（3）计算概率

本文采用AIC（Akaike Information Criterion）、BIC（Bayesian Information Criterion）法进行模型检验，AIC、BIC值越小则Copula联合分布函数的拟合程度越高，计算公式如下。

式中n为联合观测值的总数，Pei为经验频率，Pi为理论频率，k为参数个数，MSE为均方误差，RMSE为均方根误差。

式中Pe（xi，yi）为二维联合分布经验频率，Num（xj≤xi，yj≤yi）表示联合观测值小于等于值（xi，yi）的个数。

式中Pe（xi，yi，zi）为三维联合分布经验频率，Num（xj≤xi，yj≤yi，zj≤zi）表示联合观测值小于等于值（xi，yi，zi）的个数。
（4）计算重现期

4 代码片段

复现《MATLAB统计分析与应用：40个案例分析》中第七章节案例18：沪深股市日收益率的二元Copula模型（仅为示例代码，完整代码获取方式见文末）：

%******求Copula中参数的估计值******
rho_norm = copulafit('Gaussian',[U(:), V(:)])
[rho_t,nuhat,nuci] = copulafit('t',[U(:), V(:)])
linspace(x1,x2,N)
%******绘制Copula的密度函数和分布函数图******
[Udata,Vdata] = meshgrid(linspace(0,1,31));
Cpdf_norm = copulapdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);
Ccdf_norm = copulacdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);
Cpdf_t = copulapdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);
Ccdf_t = copulacdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);
% 绘制二元正态Copula的密度函数和分布函数图
figure;
surf(Udata,Vdata,reshape(Cpdf_norm,size(Udata)));
xlabel('U');
ylabel('V');
zlabel('c(u,v)');
figure;
surf(Udata,Vdata,reshape(Ccdf_norm,size(Udata)));
xlabel('U');
ylabel('V');
zlabel('C(u,v)');
% 绘制二元t-Copula的密度函数和分布函数图
figure;
surf(Udata,Vdata,reshape(Cpdf_t,size(Udata)));
xlabel('U');
ylabel('V');
zlabel('c(u,v)');
figure;
surf(Udata,Vdata,reshape(Ccdf_t,size(Udata)));
xlabel('U');
ylabel('V');
zlabel('C(u,v)');
% 绘制经验Copula分布函数图像
surf(Udata,Vdata,reshape(CopulaEmpirical,size(Udata)))
xlabel('U');
ylabel('V');
zlabel('Empirical Copula C(u,v)');

如果需要更多关于二维copula的联合重现期、同现重现期和条件重现期等的代码（MATLAB和R语言），如下图所示。详见链接（代码可直接运行，替换数据即可运行出结果）https://mbd.pub/o/bread/ZZWWlZ5t
欢迎关注公众号，获取更多科研代码和前沿论文资讯等相关内容

5 参考资料

`
[1] Li W, Jiang S, Zhao Y, Wang H, et al. A copula-based security risk evaluation and probability calculation for water-energy-food nexus[J]. Science of The Total Environment, 2023, 856: 159236.

[2] 李溦,姜珊,赵勇等.京津冀水-能源-粮食耦合系统安全评价[J].水资源保护,2023,39(05):39-48.

[3] 许怡然,鲁帆,戴雁宇等.大清河流域气象干旱时空演化和联合概率分析[J].人民黄河,2021,43(03):84-89.

[4] 宋松柏.Copulas函数及其在水文中的应用[M].科学出版社, 2012.

[5] 王丽芳.Copula分布估计算法[M].机械工业出版社, 2012.

[6] 陈璐.Copula函数理论在多变量水文分析计算中的应用研究[M].武汉大学出版社,2013.

[7] 李霞.COPULA方法及其应用[M].经济管理出版社, 2014.

[8] 谢中华.MATLAB统计分析与应用:40个案例分析[M].北京航空航天大学出版社, 2015.

[9] Chen, S., Xu, J., Li, Q., Tan, X., Nong, X., 2019. A copula-based interval-bistochastic programming method for regional water allocation under uncertainty. Agric. Water Manag. 217, 154–164. https://doi.org/10.1016/j.agwat.2019.02.008

[10] Huang, F., Ochoa, C.G., 2022. A copula incorporated cellular automata module for modeling the spatial distribution of oasis recovered by ecological water diversion: An application to the Qingtu Oasis in Shiyang River basin, China. J. Hydrol. 608, 127573. https://doi.org/10.1016/j.jhydrol.2022.127573

[11] Jin, H., Chen, X., Zhong, R., Pan, Y., Zhao, T., Liu, Z., Tu, X., 2022. Joint probability analysis of water and sediment and predicting sediment load based on copula function. Int. J. Sediment Res. 37, 639–652. https://doi.org/10.1016/j.ijsrc.2022.04.002

[12] Ren, K., Huang, S., Huang, Q., Wang, H., Leng, G., Fang, W., Li, P., 2020. Assessing the reliability, resilience and vulnerability of water supply system under multiple uncertain sources. J. Clean. Prod. 252, 119806. https://doi.org/10.1016/j.jclepro.2019.119806

[13] Ribeiro, A.F.S., Russo, A., Gouveia, C.M., Páscoa, P., 2019. Copula-based agricultural drought risk of rainfed cropping systems. Agric. Water Manag. 223, 105689. https://doi.org/10.1016/j.agwat.2019.105689

[14] Tao, Y., Wang, Y., Wang, D., Ni, L., Wu, J., 2021. A C-vine copula framework to predict daily water temperature in the Yangtze River. J. Hydrol. 598, 126430. https://doi.org/10.1016/j.jhydrol.2021.126430

[15] Zhang, Y., Guo, L., Liang, C., Zhao, L., Wang, J., Zhan, C., Jiang, S., 2022. Encounter risk analysis of crop water requirements and effective precipitation based on the copula method in the Hilly Area of Southwest China. Agric. Water Manag. 266, 107571. https://doi.org/10.1016/j.agwat.2022.107571

[16] 杨志勇,袁喆,方宏阳等.基于Copula函数的滦河流域旱涝组合事件概率特征分析[J].水利学报,2013,44(05):556-561+569.DOI:10.13243/j.cnki.slxb.2013.05.007.

[17] 杨丰顺,邵东国,顾文权等.基于Copula函数的区域需水量随机模拟[J].农业工程学报,2012,28(18):107-112.

[18] 涂新军, 陈晓宏, 赵勇, 等. 变化环境下东江流域水文干旱特征及缺水响应[J]. 水科学进展, 2016, 27(6): 810-821.

[19] 郭嘉豪, 王会肖, 赵茹欣, 等. 基于多维Copula函数的松嫩平原干旱特征分析[J]. 节水灌溉, 2020(12): 57-63+67.

[20] 黄强, 陈子燊. 基于二次重现期的多变量洪水风险评估[J]. 湖泊科学, 2015, 27(2): 352-360.

[21] 肖名忠, 张强, 陈永勤, 等. 基于三变量Copula函数的东江流域水文干旱频率分析[J]. 自然灾害学报, 2013, 22(2): 99-108.

[22] 孔祥铭. 基于Copula分析的水资源系统风险评估与管理方法研究[D]. 华北电力大学(北京), 2016.

[23] 林小敏. 基于Copula函数的灌区水资源短缺风险分析[D]. 郑州大学, 2017.

[24] 曾文颖, 徐明庆, 宋松柏, 等. 基于R-Vine copula函数的极端降水联合分布及风险识别[J]. 水资源保护: 1-12.

[25] 杨蕊, 吴时强, 高学平, 等. 基于Vine Copula函数的河湖连通水环境多因子联合风险识别研究[J]. 水利学报, 2020, 51(5): 606-616.

[26] 王燕. 水与能源耦合模拟与综合风险评估研究[D]. 中国水利水电科学研究院, 2020.
[27] 科研小技巧 | 多变量联合分布函数Copula. https://mp.weixin.qq.com/s/KmssU9cEaNE_fktThAZNMw