コピュラとは・・・

多変量確率分布を周辺分布の掛け算で求められるようにするためにデータ間の関係(相関)など諸々詰め込んで表現した関数(という理解)
https://en.wikipedia.org/wiki/Copula_(probability_theory)
多変量確率分布のパラメトリック推定に使えるらしい。実際に使ってみよう。

コピュラのモジュールをインストール

これ→ https://pypi.org/project/copulas/

In [27]:
#インストール
#!pip install copulas

Copulasに掲載されているテストコード

乱数作って正規分布でフィッティングした結果を比較する。

In [6]:
import warnings
warnings.filterwarnings('ignore')

from copulas.datasets import sample_trivariate_xyz
from copulas.multivariate import GaussianMultivariate
from copulas.visualization import compare_3d
%matplotlib inline

# Load a dataset with 3 columns that are not independent
#なんか適当なデータみたい
real_data = sample_trivariate_xyz()

# Fit a gaussian copula to the data
copula = GaussianMultivariate()
copula.fit(real_data)

# Sample synthetic data
synthetic_data = copula.sample(len(real_data))

# Plot the real and the synthetic data to compare
compare_3d(real_data, synthetic_data)

テストコードで作成した分布は正規分布っぽくなくて,うまく推定できていないように思える。元の分布が推定したい分布らしくなければうまく推定できない(あたりまえ)
だから,正規分布の乱数を作って推定してみよう。うまくいくはず

3次元の正規分布に従う乱数を生成

In [7]:
import general_module as gm #自作モジュール。乱数生成
import numpy as np
import pandas as pd
In [8]:
#次元数
dim     = 3
#学習対象の分布
#平均
mu_ar   = np.array([3,5,0])
#標準偏差
sig_ar  = np.array([2,1,3])
#相関係数(対象行列なので、作り方を考える。対角成分は自己相関なので1
rho_ar  = np.array([[1,0,0],[0,1,0],[0,0,1]])
In [28]:
#多次元分布の共分散行列を作る。これが乱数生成のタネとなる
cov_mat = np.array([])
d1=0
while d1 < dim:
    d2=0
    line=np.array([])
    #行を作る
    while d2 < dim:
        line=np.append(line, rho_ar[d1][d2]*sig_ar[d1]*sig_ar[d2])
        d2+=1
    if(d1==0):#行を作る
        cov_mat=line
    else:#行列を連結する
        cov_mat=np.vstack([cov_mat, line])

    d1+=1
In [29]:
#作った行列の確認
cov_mat
Out[29]:
array([[4., 0., 0.],
       [0., 1., 0.],
       [0., 0., 9.]])
In [12]:
#乱数シードと生成乱数の数
seed=1
num_rand=1000
In [16]:
#多次元同時正規分布を作る
import pandas as pd
multi_rand_tuple = pd.DataFrame(gm.generate_multivariate_normal_upper3d(mu_ar, cov_mat, num_rand, seed))
In [30]:
#ガウシアンコピュラでフィッティング(推定)
# Fit a gaussian copula to the data
copula = GaussianMultivariate()
copula.fit(multi_rand_tuple)

# Sample synthetic data
synthetic_data = copula.sample(len(multi_rand_tuple))

# Plot the real and the synthetic data to compare
compare_3d(multi_rand_tuple, synthetic_data)

Real DataとSynthetic Dataを重ねてみよう

In [22]:
#3Dscatterをプロット。形を見てみる
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

fig = plt.figure(figsize=(10,10))

ax = fig.add_subplot(111, projection='3d')
ax.scatter(multi_rand_tuple[0].values, multi_rand_tuple[1].values, multi_rand_tuple[2].values,alpha=0.2,color='red')
ax.scatter(synthetic_data[0].values, synthetic_data[1].values, synthetic_data[2].values,alpha=0.1,color='blue')

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

あってるっぽいね(定量評価は今回はスルーします)