Up

LKJ相関行列確率分布シミュレーション

 

多変量正規分布は基本的な確率分布の１つであり、共分散行列をパラメータとしてもつ。共分散行列は標準偏差と相関行列により表すことができるが、相関行列の確率分布としてLKJ相関行列確率分布（LKJ correlation distribution）がある。Lambert (2018) “A student’s guide to Bayesian Statistics” は、相関行列の確率分布としてWishart分布は歴史的な事情で用いられるもので、現代ではLKJ相関行列確率分布を用いることを勧めている。

LKJ相関行列確率分布はパラメータを１つもち、確率分布は次式で与えられる。

このとき、のサイズを、の要素をとするとき、

である。

LKJ相関行列確率分布に従う乱数を発生するStanスクリプトを以下のように用意した。

 

    generated quantities {

        corr_matrix[3] corr;

        corr = lkj_corr_rng(3, %s);   

    }

 

このStanスクリプトを用いるPythonスクリプトはリスト１のように用意した。なお、本ウェブサイトのプログラムは、圧縮ファイルLKJSimuFiles.zipにまとめた。

モデル宣言部で乱数生成を宣言する場合は、以下のようになる。

 

    parameters {

        corr_matrix[3] corr;

    }

    model {

        corr ~ lkj_corr(%s);

    }

 

モデル宣言部で乱数生成のスクリプトを用意した場合のPythonスクリプトはリスト２に示した。

リスト１あるいはリスト２のプログラムを実行すると、の値の入力を求められので、適当な値（正の値）を以下のように入力する。

 

(py39) PS D:\LKJSimuFiles> python .\LKJCorrSimuGenQuant.py

eta = 1.0

 

「eta =」の右側に適当な値（正の値）を入力して、Enterキーを押すと計算が始まる。

乱数の生成が終了すると、まず相関係数のペアごとの散布図が表示される（図１.1）。

図1.1

 

図1.1のフォームを閉じると、各相関係数のヒストグラムが表示される（図1.2）。

図1.2

 

各相関係数は、式（1.2）のベータ分布に従う。

式（1.1）の右辺は、相関行列の行列式の乗であるが、のときは０乗であり定数１である。

図1.2のフォームを閉じると、相関行列の行列式のヒストグラムが表示される（図1.3）。

図1.3

 

相関行列の行列式は、一様分布ではないことが分かる。式（1.1）の意味するところは、一様分布しているのは相関行列であって、その行列式ではないということである。このことは、次の例を見れば明らかである。

 

分布が原点を除いた単位円内で、確率密度が原点からの距離の乗に比例している次式の分布を考える。

 

 

式（2.1）に従うサンプルは、次のMetropolis algorithmで生成することができる。

 

def mcmc_sampling(eta = 1.0, n = 1000):

    """

            Metropolis algorithm

    """

    def p_dist(x, y):

        v = (x**2 + y**2) ** 0.5

        if v <= 0.0 or v > 1.0:

            p = 0.0

        else:

            p = v ** (eta - 1.0)

        return p

 

    def next_xy(x0, y0):

        x1 = np.random.normal(x0, 0.3)

        y1 = np.random.normal(y0, 0.3)

        t_p = p_dist(x1, y1)

        if t_p <= 0.0:

            x1 = x0

            y1 = y0

        else:

            a = t_p / p_dist(x0, y0)

            if np.random.uniform() > a:

                x1 = x0

                y1 = y0

        return x1, y1

 

    x0 = 0.5

    y0 = 0.5

    for i in range(100):

        x0, y0 = next_xy(x0, y0)

    Xs = []

    Ys = []

    xt = x0

    yt = y0

    for i in range(n):

        xt, yt = next_xy(xt, yt)

        Xs.append(xt)

        Ys.append(yt)

 

    return Xs, Ys

 

この式（2.1）に従うサンプルを得る関数によりサンプリングを行うPythonスクリプトをリスト３に示す。スクリプトファイルは、圧縮ファイルLKJSimuFiles.zipにまとめている。

リスト３のスクリプトを実行するとの値の入力を求めてくる。

 

(py39) PS D:\LKJSimuFiles> python CircularDist.py

eta = 1

 

上の例では、「eta = 1」と値１を入力している。入力後、Enterキーを押すと計算が始まる。

計算が終了すると図2.1が表示される。

図2.1

 

左上の散布図を見ると、単位円内に一様に分布していることが分かる。式（2.1）は、のときは距離に関わらず一様分布であることを意味しているが、距離のヒストグラムは右肩上がりのグラフである。これは、円周の長さが、半径に比例していることによる。このことは、LKJ分布（1.1）がのとき行列式に関して定数であるが、行列式のヒストグラム（図1.3）は一様分布でないことと対応している。図2.1の下段に周辺分布のヒストグラムが示されているが、これも一様分布ではない。

 

の場合のグラフは図2.2のようになる。点が原点に集まっていることが分かる。

図2.2

 

とすると、図2.3のグラフを得る。分布が原点から遠ざかり周辺に集まっていることが分かる。

図2.3

 

PythonとPyStanの説明は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。

 

 

リスト１　LKJ相関行列確率分布生成のPythonスクリプト（Stanのサンプル生成部を用いる場合；ファイル名LKJCorrSimuGenQuant.py）

 

import pystan

import matplotlib.pyplot as plt

 

StanCode = """

    generated quantities {

        corr_matrix[3] corr;

        corr = lkj_corr_rng(3, %s);   

    }

"""

eta = input('eta = ')

sm = pystan.StanModel(model_code = (StanCode % eta))

 

fit = sm.sampling(algorithm="Fixed_param", n_jobs = 1) 

print(fit)

 

r12 = []

r13 = []

r23 = []

for r in fit['corr']:

    r12.append(r[0,1])

    r13.append(r[0,2])

    r23.append(r[1,2])

 

plt.subplots(figsize = (7,7))

plt.subplot(2,2,1)

plt.plot(r12, r13, 'o', alpha = 0.3)

plt.xlabel('r12\n', fontsize = 16)

plt.ylabel('\nr13', fontsize = 16)

plt.title('$\eta$ = {}'.format(eta))

 

plt.subplot(2,2,2)

plt.plot(r23, r13, 'o', alpha = 0.3)

plt.xlabel('r23\n', fontsize = 16)

plt.ylabel('\nr13', fontsize = 16)

plt.title('$\eta$ = {}'.format(eta))

 

plt.subplot(2,2,4)

plt.plot(r23, r12, 'o', alpha = 0.3)

plt.xlabel('r23\n', fontsize = 16)

plt.ylabel('\nr12', fontsize = 16)

plt.title('$\eta$ = {}'.format(eta))

 

plt.tight_layout()

plt.show()

 

 

plt.subplots(figsize = (5,9))

plt.subplot(3,1,1)

plt.hist(r12)

plt.title(' ', fontsize = 18)

plt.xlabel('r12', fontsize = 18)

plt.title('$\eta$ = {}'.format(eta))

 

plt.subplot(3,1,2)

plt.hist(r13)

plt.title(' ', fontsize = 18)

plt.xlabel('r13', fontsize = 18)

plt.title('$\eta$ = {}'.format(eta))

 

plt.subplot(3,1,3)

plt.hist(r23)

plt.title(' ', fontsize = 18)

plt.xlabel('r23', fontsize = 18)

plt.title('$\eta$ = {}'.format(eta))

 

plt.tight_layout()

plt.show()

 

import numpy as np

 

det_Rs = []

for R in fit['corr']:

    det_Rs.append(np.linalg.det(R))

   

plt.hist(det_Rs)

plt.xlabel('Det of Corr Matrix')

plt.title('$\eta$ = {}'.format(eta))

plt.show()

 

 

リスト２　モデル宣言部で乱数を生成するStanスクリプトを利用したPythonスクリプト（ファイル名LKJCorrSimuModel.py）

 

import pystan

import matplotlib.pyplot as plt

 

StanCode = """

    parameters {

        corr_matrix[3] corr;

    }

    model {

        corr ~ lkj_corr(%s);

    }

"""

eta = input('eta = ')

 

sm = pystan.StanModel(model_code = (StanCode % eta))

fit = sm.sampling(n_jobs = 1)

print(fit)

 

r12 = []

r13 = []

r23 = []

for r in fit['corr']:

    r12.append(r[0,1])

    r13.append(r[0,2])

    r23.append(r[1,2])

 

plt.subplots(figsize = (7,7))

plt.subplot(2,2,1)

plt.plot(r12, r13, 'o', alpha = 0.3)

plt.xlabel('r12\n', fontsize = 16)

plt.ylabel('\nr13', fontsize = 16)

plt.title('$\eta$ = {}'.format(eta))

 

plt.subplot(2,2,2)

plt.plot(r23, r13, 'o', alpha = 0.3)

plt.xlabel('r23\n', fontsize = 16)

plt.ylabel('\nr13', fontsize = 16)

plt.title('$\eta$ = {}'.format(eta))

 

plt.subplot(2,2,4)

plt.plot(r23, r12, 'o', alpha = 0.3)

plt.xlabel('r23\n', fontsize = 16)

plt.ylabel('\nr12', fontsize = 16)

plt.title('$\eta$ = {}'.format(eta))

 

plt.tight_layout()

plt.show()

 

 

plt.subplots(figsize = (5,9))

plt.subplot(3,1,1)

plt.hist(r12)

plt.title(' ', fontsize = 18)

plt.xlabel('r12', fontsize = 18)

plt.title('$\eta$ = {}'.format(eta))

 

plt.subplot(3,1,2)

plt.hist(r13)

plt.title(' ', fontsize = 18)

plt.xlabel('r13', fontsize = 18)

plt.title('$\eta$ = {}'.format(eta))

 

plt.subplot(3,1,3)

plt.hist(r23)

plt.title(' ', fontsize = 18)

plt.xlabel('r23', fontsize = 18)

plt.title('$\eta$ = {}'.format(eta))

 

plt.tight_layout()

plt.show()

 

import numpy as np

 

det_Rs = []

for R in fit['corr']:

    det_Rs.append(np.linalg.det(R))

   

plt.hist(det_Rs)

plt.xlabel('Det of Corr Matrix')

plt.title('$\eta$ = {}'.format(eta))

plt.show()

 

 

リスト３　環状分布の例（ファイル名CircularDist.py）

 

import numpy as np

import matplotlib.pyplot as plt

 

 

def mcmc_sampling(eta = 1.0, n = 1000):

    """

            Metropolis algorithm

    """

    def p_dist(x, y):

        v = (x**2 + y**2) ** 0.5

        if v <= 0.0 or v > 1.0:

            p = 0.0

        else:

            p = v ** (eta - 1.0)

        return p

 

    def next_xy(x0, y0):

        x1 = np.random.normal(x0, 0.3)

        y1 = np.random.normal(y0, 0.3)

        t_p = p_dist(x1, y1)

        if t_p <= 0.0:

            x1 = x0

            y1 = y0

        else:

            a = t_p / p_dist(x0, y0)

            if np.random.uniform() > a:

                x1 = x0

                y1 = y0

        return x1, y1

 

    x0 = 0.5

    y0 = 0.5

    for i in range(100):

        x0, y0 = next_xy(x0, y0)

    Xs = []

    Ys = []

    xt = x0

    yt = y0

    for i in range(n):

        xt, yt = next_xy(xt, yt)

        Xs.append(xt)

        Ys.append(yt)

 

    return Xs, Ys

 

eta = float(input('eta = '))

 

Xs, Ys = mcmc_sampling(eta = eta, n = 10_000)

 

plt.subplots(figsize = (8, 8))

plt.subplot(2,2,1)

plt.plot(Xs, Ys, 'o', alpha = 0.3)

plt.xlabel('X', fontsize = 14)

plt.ylabel('Y', fontsize = 14)

plt.title('$\eta$ = {0:.2f}'.format(eta), fontsize = 16)

 

Ds = []

for i in range(len(Xs)):

    Ds.append(((Xs[i]**2) + (Ys[i]**2))**0.5)

plt.subplot(2,2,2)

plt.hist(Ds)

plt.xlabel('Distance', fontsize = 14)

plt.title('Histogram of Distances', fontsize = 16)

 

plt.subplot(2,2,3)

plt.hist(Xs)

plt.xlabel('X', fontsize = 14)

plt.title('\nHistogram of X', fontsize = 16)

 

plt.subplot(2,2,4)

plt.hist(Ys)

plt.xlabel('Y', fontsize = 14)

plt.title('\nHistogram of Y', fontsize = 16)

 

plt.tight_layout()

plt.show()

 

 

 

Up