それぞれが任意の分布に従う，相関を持つ２確率変数の（近似的）作り方メモ

要点

• 標準正規分布に従う確率変数X，χ2分布に従う確率変数Yを，任意の相関係数ρで作成する方法を探していました．
• X, Yともに標準正規分布に従うならば，以下の方法で作成することができます．
  • numpy.random.multivariate_normal関数を使用する．
    • 多次元標準正規分布ならば，分散共分散行列が相関行列と等しくなるので, covmtx=[[1,ρ],[ρ,1]]
• 標準正規分布でなくても，平均がゼロで分散が等しい分布ならば，horiemさん解説のこの方法で生成できます．
• X~標準正規分布とY~χ2分布はこの方法の条件を満たさないので、別の方法で生成する必要があります．
• 一時しのぎの荒業として，多次元標準正規分布で相関を持たせて2確率変数を生成し，その片方を別の確率分布に従うように変換する方法をメモしておきます．
  • Y -> 標準正規分布の累積分布cdf(Y) -> 任意の分布のパーセント点関数ppf(cdf(Y)) -> Z
  • なんのこっちゃ、というときはこちらを参照のこと．
  • 変換後の相関係数の値を厳密に指定できないという欠点があるので、要注意．

コード

sample.py

# import numpy as np
# import scipy.stats as spstats

# generate X,Y ~ multinormal(mu = [0,0], cov = [[1,ρ],[ρ,1]])
rho_norm = 0.8 # correlation coeff for multinorm
mu = [0, 0] # mean of X, Y
cov = [
         [1, rho_norm], 
         [rho_norm, 1]
      ] # cov matrix

vals_norm = np.random.multivariate_normal(mu, cov, 100000)
x_norm = vals_norm[:,0]
y_norm = vals_norm[:,1]
# np.corrcoef(x_norm, y_norm) gives a rho value around rho_norm

# convert Y to Z ~ chi2(k)
k = 3 # parameter of chi2 dist
z_chi2 = spstats.chi2.ppf(spstats.norm.cdf(y_norm, loc = 0, scale = 1), df = k)

# x_norm ~ norm(mu = 0, var = 1) and z_chi2 ~ chi2(k = 3)
# np.corrcoef(x_norm, z_chi2) gives a rho value a bit smaller than rho_norm

結果

相関係数0.8で生成したX, Yの散布図と周辺分布はこんな感じになります．

で，上記の変換を噛ませたZについて，Xとの散布図と周辺分布はこんな感じになります．

Z(縦軸)について，周辺分布がk=3のχ2分布に変換されていることがわかります．
このときの相関係数（相関行列が求まるので，[0,1]成分）は

正の相関は保たれていますが，多次元標準正規分布でX-Y間に指定した相関係数(0.8)より小さくなっちゃってますね．# -> 2021/12/15 コメントを受け末尾註追記

今回はYをχ2分布へ変換しましたが，任意の分布へ変換することが可能です．また，Xも同様に変換可能です．

注意点

この方法では変換後の相関係数を厳密に指定することができません．
もっと正攻法は方法があるように思いますが，ちょっと見つけることができなかったので，ここにメモしておきます．
誰かもっとやりやすい方法教えてください．

11/26追記：matlabの本家解説にこんなのありました．ほぼ同じアプローチですが，より詳細に記述されています．FMI(For My Info)

2021/12/15追記：非線形な関係に変形しているので、ピアソンの相関係数が小さくなるのは当たり前でしたね。コメントで指摘いただいたように、スピアマンの順位相関係数などで評価をするなら、相関係数は大きく変化していないはずです。
関連1:Pythonで相関がある任意の確率分布に従う乱数を生成する
関連2:Pythonでモンテカルロ法を用いた建設コストリスクの定量化（多変量解析） <- 多変量で同様の変換を実施