python中Copula在多元联合分布建模可视化2实例合集|附数据代码
相关视频
这样，变量的随机性和耦合性就被分离开来，其中随机变量各自的随机性由边缘分布描述，而随机变量之间的耦合特性则由Copula函数描述。
本文旨在通过一系列实例，展示如何在Python中使用Copula进行多元联合分布建模和可视化。我们将从简单的二元Copula模型开始，逐步过渡到更复杂的多元模型，并介绍如何使用不同的Copula类型和参数来适应不同的数据特性。
1.Copula在多元联合分布建模
Copula函数在金融风险管理、精算学和统计推断等领域有广泛应用。它几乎包含了随机变量所有的相依信息，因此对于分析变量之间的相关关系非常有用，尤其是在传统的线性相关系数可能无法准确度量相关关系的情况下。
具体来说，Copula函数是一个从[0,1]^n到[0,1]的映射，用于链接n个随机变量的边缘累积分布函数。它用于描述多元随机变量之间的依赖关系，这些关系可以是正相关、负相关或无相关。
在建模系统时，经常会遇到涉及多个参数的情况。这些参数中的每一个都可以用给定的概率密度函数（PDF）来描述。如果想要生成一组新的参数值，就需要从这些分布（也称为边缘分布）中进行抽样。主要有两种情况：?(i)? PDF是独立的；?(ii)??存在依赖关系。一种建模依赖关系的方法是使用copula。
Kaizong Ye
拓端分析师
从copula中抽样?
让我们用一个双变量示例来说明，并假设首先我们有一个先验知识，并知道如何对两个变量之间的依赖关系进行建模。
在这种情况下，我们使用Gumbel copula并固定其超参数theta=2。我们可以可视化其二维PDF。
_ = copulot_pdf() # 可视化
并且我们可以从该PDF中进行抽样。
现在让我们再回到这两个变量上来。在这种情况下，我们考虑它们是服从伽马分布和正态分布的。如果它们彼此独立，我们可以单独从每个PDF中进行抽样。这里我们使用一个方便的类来执行相同的操作。
左右滑动查看更多
01
02
03
04
可重复性?
_ = h.set_labels("X1", "X2", fontsize=16)
现在，我们已经使用copula表达了变量之间的依赖关系，我们可以使用同一个方便的类从这个copula中抽样得到一组新的观测值。
# 使用一个已初始化的Generator对象
h = snjoind="scatter")
这里有两点需要注意。(i)就像独立情况一样，边缘分布正确地显示了伽马分布和正态分布；(ii)两个变量之间的依赖关系可见。
估计copula参数?
现在，假设我们已经有了实验数据，并且知道可以使用Gumbel copula来表达依赖关系。但是我们不知道copula的超参数值是多少。在这种情况下，我们可以估计这个值。
我们将使用刚才生成的样本，因为我们已经知道应该得到的超参数值：theta=2。
fit_carm(sample)
print(theta)
我们可以看到，估计的超参数值接近之前设定的值。
2.python中的copula：Frank、Clayton和Gumbel copula模型估计与可视化
这篇文章中即将出现的大部分内容都会用Jupyter Notebooks来构建。
软件
scikit-learn或scipy中没有明确的copula包的实现。
2D数据的Frank、Clayton和Gumbel copula
测试
第一个样本（x）是从一个β分布中产生的，（y）是从一个对数正态中产生的。β分布的支持度是有限的，而对数正态的右侧支持度是无穷大的。对数的一个有趣的属性。两个边际都被转换到了单位范围。
我们对样本x和y拟合了三个族（Frank, Clayton, Gumbel）的copulas，然后从拟合的copulas中提取了一些样本，并将采样输出与原始样本绘制在一起，以观察它们之间的比较。
#等同于ppf，但直接从数据中构建
sortedvar=np.sort(var)
#绘制
for index,family in enumerate(['Frank', 'clayton', 'gumbel']):
#获得伪观测值
u,v = copula_f.generate_uv(howmany)
#画出伪观测值
axs[index][0].scatter(u,v,marker='o',alpha=0.7)
plt.show()
#总样本与伪观测值的对比
sz=300
loc=0.0 #对大多数分布来说是需要的
sc=0.5
y=lognorm.rvs(sc,loc=loc, size=sz)
独立（不相关）数据
我们将从β分布中抽取（x）的样本，从对数正态中抽取（y）的样本。这些样本是伪独立的（我们知道，如果你用计算机来抽取样本，就不会有真正的独立，但好在是合理的独立）。
#不相关的数据：一个β值（x）和一个对数正态（y）。
a= 0.45#2. #alpha
b=0.25#5. #beta
#画出不相关的x和y
plt.plot(t, beta.pdf(t,a,b), lw=5, alpha=0.6, label='x:beta')
#绘制由不相关的x和y建立的共线性图
title='来自不相关数据的共线性 x: beta, alpha {} beta {}, y: lognormal, mu {}, sigma dPlot(title,x,y,pseudoobs)
相依性（相关）数据
自变量将是一个对数正态（y），变量（x）取决于（y），关系如下。初始值为1（独立）。然后，对于每一个点i, 如果
, 那么?, 其中c是从1的分数列表中统一选择的，否则,?.
#相关数据：一个对数正态（y）。
#画出相关数据
np.linspace(0, lognorm.ppf(0.99, sc), sz)
plt.plot(t, gkxx.pdf(t), lw=5, alpha=0.6,
拟合copula参数
没有内置的方法来计算archimedean copulas的参数，也没有椭圆elliptic copulas的方法。但是可以自己实现。选择将一些参数拟合到一个scipy分布上，然后在一些样本上使用该函数的CDF方法，或者用一个经验CDF工作。这两种方法在笔记本中都有实现。
因此，你必须自己写代码来为archimedean获取参数，将变量转化为统一的边际分布，并对copula进行实际操作。它是相当灵活的。
#用于拟合copula参数的方法
# === Frank参数拟合
"""
对这个函数的优化将给出参数
"""
#一阶debye函数的积分值 int_debye = lambda t: t/(npexp(t)-1.)
debye = lambda alphaquad(int_debye ,
alpha
)[0]/alpha
diff = (1.-kTau)/4.0-(debye(-alpha)-1.)/alpha
#================
#clayton 参数方法
def Clayton(kTau):
try:
return 2.*kTau/(1.-kTau)
#Gumbel参数方法
def Gumbel(kTau):
try:
return 1./(1.-kTau)
#================
#copula生成
#得到协方差矩阵P
#x1=norm.ppf(x,loc=0,scale=1)
#y1=norm.ppf(y,loc=0,scale=1)
#return norm.cdf((x1,y1),loc=0,scale=P)
#================
#copula绘图
fig = pylab.figure()
ax = Axes3D(fig)
ax.text2D(0.05, 0.95, label, transform=ax.transAxes)
ax.set_xlabel('X: {}'.format(xlabel))
ax.set_ylabel('Y: {}'.format(ylabel))
#sample是一个来自U,V的索引列表。这样，我们就不会绘制整个copula曲线。
if plot:
print "绘制copula {}的样本".format(copulaName)
returnable[copulaName]=copulapoints
if plot:
zeFigure=plot3d(U[样本],V[样本],copulapoints[样本], label=copulaName,
生成一些输入数据
在这个例子中，我们使用的是与之前相同的分布，探索copula 。如果你想把这段代码改编成你自己的真实数据，。
t = np.linspace(0, lognorm.ppf(0.99, sc), sz)
#从一些df中抽取一些样本
X=beta.rvs(a,b,size=sz)
Y=lognorm.rvs(sc,size=sz)
#通过对样本中的数值应用CDF来实现边缘分布
U=beta.cdf(X,a,b)
V=lognorm.cdf(Y,sc)
#画出它们直观地检查独立性
plt.scatter(U,V,marker='o',alpha=0.7)
plt.show()
可视化Copulas
没有直接的构造函数用于高斯或tCopulas，可以为椭圆Copulas(EllipticCopulas)建立一个更通用的函数。
Frechét-H?ffding边界可视化
根据定理，我们将copula画在一起，得到了Frechét-H?ffding边界。
#建立边界为copula的区域
plot_trisurf(U[样本],V[样本],copula['min'][样本],
c='red') #上限
plot_trisurf(U[样本],V[样本],copula['max'][样本],
c='green') #下限
- 发表于:
- 原文链接：https://page.om.qq.com/page/OmhCmaQLZeCWANpyMkUdR4XQ0
- 腾讯「腾讯云开发者社区」是腾讯内容开放平台帐号（企鹅号）传播渠道之一，根据《腾讯内容开放平台服务协议》转载发布内容。
- 如有侵权，请联系 cloudcommunity@tencent.com 删除。
相关快讯
扫码
添加站长 进交流群
领取专属 10元无门槛券
私享最新 技术干货
腾讯云开发者
