Up

2変量正規分布によるベイズ分析

StanPython

 

2つの変量の対でデータが得られた場合を考える。20人の学生の国語と英語の点数のデータが図1のようであったとする。

図1 Excelデータ(TwoSubj.xlsx

 

このデータに対して2変量正規分布モデルによるベイズ分析を以下のように行ってみた。

2変量が次の2変量正規分布に従うとする。

上のモデルに対するStanスクリプトを本ウェブサイト後半に示すリスト1のように作成した。スクリプトファイルなどは、圧縮ファイルfilesTwoVars_nfd.zipにまとめているので、ダウンロード解凍すれば自由に利用できる。

リスト1のStanスクリプトによってベイズ分析を行うPythonスクリプトを本ウェブサイト後半に示すリスト2のように作成した。

リスト2のスクリプトを実行すると、まず、入力データファイル名の設定が求められる(図2)。

図2

 

図2では、図1のデータファイル名が設定されている。

入力データファイルはExcel形式のファイルを図1の形式で用意する。すなわち、1行目に変数名を置き、データは2行目以降に書く。1列目はデータの通し番号であるが、文字列でもよい。2列目と3列目にそれぞれの変数のデータ値を書く。

入力データファイル名を設定して、「Enter」キーを押すと、次に画像保存用ファイル名の設定が求められる(図2)。事後分布のグラフを保存するファイルの名前であるが、拡張子が「.png」の任意のファイル名でよい。複数の画像ファイルが保存されるが、設定ファイル名の「.png」の部分を「_mu1.png」などと変更した名前のファイルで保存される。図2では、「MyFigure.png」と設定されている。

画像保存ファイル名を設定して「Enter」キーを押すと、データが読み込まれ、処理が始まる。Stanスクリプトのコンパイルとビルドに多少時間が掛かるが、その後、MCMCサンプリングが始まる。サンプリング終了後、まず、の事後分布が表示される(図3)。

図3

 

図3のウィンドウを閉じると、の事後分布が表示される(図4)。

図4

 

図4のウィンドウを閉じると、の事後分布が表示される(図5)。

図5

 

図5のウィンドウを閉じると、の事後分布が表示される(図6)。

図6

 

図6のウィンドウを閉じると、相関係数の事後分布が表示される(図7)。

図7

 

図7のウィンドウを閉じると、プログラムの実行終了である(図8)。

図8

 

 

 

リスト1 2変量正規分布モデルのStanスクリプト(TwoVarsNormal.stan

 

data {

    int N;

    vector[2] X[N];

}

parameters {

    real mu1;

    real mu2;

    real<lower = 0.0> sgm1;

    real<lower = 0.0> sgm2;

    real<lower = -1.0, upper = 1.0> rho;

}

transformed parameters {

    vector[2] Mu;

    cov_matrix[2] Sgm;

    Mu[1] = mu1;

    Mu[2] = mu2;

    Sgm[1][1] = sgm1 * sgm1;

    Sgm[1][2] = rho * sgm1 * sgm2;

    Sgm[2][1] = Sgm[1][2];   // rho * sgm2 * sgm1;

    Sgm[2][2] = sgm2 * sgm2;

}

model {

    mu1 ~ normal(0.0, 10000.0);

    mu2 ~ normal(0.0, 10000.0);

    sgm1 ~ exponential(0.0001);

    sgm2 ~ exponential(0.0001);

    rho ~ uniform(-1.0, 1.0);

    for (i in 1:N)

        X[i] ~ multi_normal(Mu, Sgm);

}

 

 

リスト2 リスト1Stanスクリプトを用いるPythonスクリプト(MainTwoVars.py

 

import pandas as pd

import numpy as np

import pystan

import matplotlib.pyplot as plt

import seaborn as sb

 

finnm = input('入力データファイル(*.xlsx= ')

foutnm = input('画像保存ファイル(*.png= ')

#

#      Excelファイルの読み込み

#   モジュール xlrd がインストール済みであること

#

xlsx_fl = pd.ExcelFile(finnm)                       #   ExcelファイルからExcelFile型オブジェクトを生成

data_xlsx = pd.read_excel(xlsx_fl,  header = None)  #   ExcelFile型オブジェクトからDataFrame型オブジェクトを生成

data = data_xlsx.values                             #   Dataframe型オブジェクトからnumpyの配列オブジェクトを生成

print('\n', data)

 

VarName1 = data[0, 1]

VarName2 = data[0, 2]

print(VarName1, VarName2)

X = data[1:, 1:3]

X = [[v[0], v[1]] for v in X]

print(X)

   

DataStan = {'N': len(X), 'X': X}

 

sm = pystan.StanModel(file = 'TwoVarsNormal.stan')

fit = sm.sampling(data = DataStan, n_jobs = 1)

print(fit)

 

SavedFiles = []

 

#       日本語表記の設定

plt.rcParams['font.family'] = ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']

 

mu1_025, mu1_50, mu1_975 = np.percentile(fit['mu1'], [2.5, 50, 97.5])

plt.title(('{0:}の平均値の事後分布'.format(VarName1)) + \

          '\nMed = {0:.2f},  95%CI = [{1:.2f}, {2:.2f}]'.format(

              mu1_50, mu1_025, mu1_975), fontsize = 20)

sb.kdeplot(fit['mu1'])

plt.xlabel(VarName1 + 'の平均値', fontsize = 16)

plt.yticks([])

plt.tight_layout()

plt.savefig(foutnm[:-4] + '_mu.png')

SavedFiles.append(foutnm[:-4] + '_mu1.png')

plt.show()

 

mu2_025, mu2_50, mu2_975 = np.percentile(fit['mu2'], [2.5, 50, 97.5])

plt.title(('{0:}の平均値の事後分布'.format(VarName2)) + \

          '\nMed = {0:.2f},  95%CI = [{1:.2f}, {2:.2f}]'.format(

              mu2_50, mu2_025, mu2_975), fontsize = 20)

sb.kdeplot(fit['mu2'])

plt.xlabel(VarName2 + 'の平均値', fontsize = 16)

plt.yticks([])

plt.tight_layout()

plt.savefig(foutnm[:-4] + '_mu2.png')

SavedFiles.append(foutnm[:-4] + '_mu2.png')

plt.show()

 

sgm1_025, sgm1_50, sgm1_975 = np.percentile(fit['sgm1'], [2.5, 50, 97.5])

plt.title(('{0:}の標準偏差の事後分布'.format(VarName1)) + \

          '\nMed = {0:.2f},  95%CI = [{1:.2f}, {2:.2f}]'.format(

              sgm1_50, sgm1_025, sgm1_975), fontsize = 20)

sb.kdeplot(fit['sgm1'])

plt.xlabel(VarName1 + 'の標準偏差', fontsize = 16)

plt.yticks([])

plt.tight_layout()

plt.savefig(foutnm[:-4] + '_sgm1.png')

SavedFiles.append(foutnm[:-4] + '_sgm1.png')

plt.show()

 

sgm2_025, sgm2_50, sgm2_975 = np.percentile(fit['sgm2'], [2.5, 50, 97.5])

plt.title(('{0:}の標準偏差の事後分布'.format(VarName2)) + \

          '\nMed = {0:.2f},  95%CI = [{1:.2f}, {2:.2f}]'.format(

              sgm2_50, sgm2_025, sgm2_975), fontsize = 20)

sb.kdeplot(fit['sgm2'])

plt.xlabel(VarName2 + 'の標準偏差', fontsize = 16)

plt.yticks([])

plt.tight_layout()

plt.savefig(foutnm[:-4] + '_sgm2.png')

SavedFiles.append(foutnm[:-4] + '_sgm2.png')

plt.show()

 

rho_025, rho_50, rho_975 = np.percentile(fit['rho'], [2.5, 50, 97.5])

plt.title('相関係数の事後分布' + \

          '\nMed = {0:.2f},  95%CI = [{1:.2f}, {2:.2f}]'.format(

              rho_50, rho_025, rho_975), fontsize = 20)

sb.kdeplot(fit['rho'])

plt.xlabel('相関係数', fontsize = 16)

plt.yticks([])

plt.tight_layout()

plt.savefig(foutnm[:-4] + '_rho.png')

SavedFiles.append(foutnm[:-4] + '_rho.png')

plt.show()

 

print('\n保存ファイル:')

for nm in SavedFiles:

    print('    ', nm)

print()

 

 

Up