Copula - 多變量聯合分佈¶
[1]:
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from scipy import stats
sns.set_style("darkgrid")
sns.mpl.rc("figure", figsize=(8, 8))
[2]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
return false;
}
在建模系統時，經常會遇到涉及多個參數的情況。每個參數都可以使用給定的機率密度函數 (PDF) 來描述。如果我們想要能夠生成一組新的參數值，我們需要能夠從這些分佈（也稱為邊緣分佈）中取樣。主要有兩種情況：(i) PDF 是獨立的；(ii) 存在依賴關係。一種對依賴關係建模的方法是使用 copula。
從 copula 中取樣¶
讓我們使用一個雙變數範例，並首先假設我們有一個先驗，並且知道如何對兩個變數之間的依賴關係進行建模。
在這種情況下，我們使用 Gumbel copula 並將其超參數 theta=2 固定。我們可以將其二維 PDF 可視化。
[3]:
from statsmodels.distributions.copula.api import (
CopulaDistribution, GumbelCopula, IndependenceCopula)
copula = GumbelCopula(theta=2)
_ = copula.plot_pdf() # returns a matplotlib figure
我們可以對 PDF 進行取樣。
[4]:
sample = copula.rvs(10000)
h = sns.jointplot(x=sample[:, 0], y=sample[:, 1], kind="hex")
_ = h.set_axis_labels("X1", "X2", fontsize=16)
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/statsmodels/tools/rng_qrng.py:54: FutureWarning: Passing `None` as the seed currently return the NumPy singleton RandomState
(np.random.mtrand._rand). After release 0.13 this will change to using the
default generator provided by NumPy (np.random.default_rng()). If you need
reproducible draws, you should pass a seeded np.random.Generator, e.g.,
import numpy as np
seed = 32839283923801
rng = np.random.default_rng(seed)"
warnings.warn(_future_warn, FutureWarning)
讓我們回到我們的兩個變數。在這種情況下，我們認為它們分別呈 gamma 分佈和常態分佈。如果它們彼此獨立，我們可以分別從每個 PDF 中取樣。在這裡，我們使用一個方便的類別來執行相同的操作。
可重複性¶
從 copula 中生成可重複的隨機值需要明確設置 seed 引數。seed 接受已初始化的 NumPy Generator 或 RandomState，或任何可接受 np.random.default_rng 的引數，例如，整數或整數序列。此範例使用整數。
直接在 np.random 分佈中公開的單例 RandomState 不會使用，且設定 np.random.seed 對生成的值沒有影響。
[5]:
marginals = [stats.gamma(2), stats.norm]
joint_dist = CopulaDistribution(copula=IndependenceCopula(), marginals=marginals)
sample = joint_dist.rvs(512, random_state=20210801)
h = sns.jointplot(x=sample[:, 0], y=sample[:, 1], kind="scatter")
_ = h.set_axis_labels("X1", "X2", fontsize=16)
現在，在上面我們使用 copula 表示了變數之間的依賴關係，我們可以使用此 copula 來取樣一組新的觀測值，並使用相同的方便類別。
[6]:
joint_dist = CopulaDistribution(copula, marginals)
# Use an initialized Generator object
rng = np.random.default_rng([2, 0, 2, 1, 0, 8, 0, 1])
sample = joint_dist.rvs(512, random_state=rng)
h = sns.jointplot(x=sample[:, 0], y=sample[:, 1], kind="scatter")
_ = h.set_axis_labels("X1", "X2", fontsize=16)
這裡有兩點需要注意。(i) 與獨立情況一樣，邊緣分佈正確地顯示了 gamma 分佈和常態分佈；(ii) 兩個變數之間存在可見的依賴關係。
估計 copula 參數¶
現在，假設我們已經有了實驗數據，並且知道可以使用 Gumbel copula 表示依賴關係。但是我們不知道 copula 的超參數值。在這種情況下，我們可以估計該值。
我們將使用剛剛生成的樣本，因為我們已經知道應該得到的超參數值：theta=2。
[7]:
copula = GumbelCopula()
theta = copula.fit_corr_param(sample)
print(theta)
2.049379621506455
我們可以發現，估計的超參數值接近先前設定的值。