Up

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

Stanの使用

 

1変量のデータ、例えば20人の学生の国語の点数が図1のように与えられているとする。

図1 Excelデータファイル(OneVar.xlsx

 

図1のデータは、1行目に変数名が書かれている。1列目はデータの識別用で、図1では通し番号が設定されているが文字列でもよい。データ値は2列目に書かれている。

このデータを正規分布モデル

によってベイズ分析するためのStanスクリプトを本ウェブサイト後半に示すリスのように用意した。スクリプトファイルなどは圧縮ファイルfilesbysonevar.zipにまとめたので、ダウンロード解凍すれば自由に用いることができる。

リスト1のStanスクリプトを利用して図1のようなデータをベイズ分析するPythonスクリプトとして、本ウェブサイト後半に示すリスを用意した。

リスト2のスクリプトを実行すると、まず入力データ名の設定を求めるダイアログウィンドウが表示される(図2)。

図2

 

図2では、図1のExcelデータファイル名が設定されている。入力データファイルは図1の形式のExcelファイル(互換ソフトで作成されたものでもよい)を指定する。

入力データファイル名を設定して「開く」ボタンをクリックすると、次に画像ファイルの保存先ファイル名の設定を求めるダイアログウィンドウが表示される(図3)。

図3

 

保存先ファイル名はファイル拡張子が「.png」であることを除いて任意である。図3では「MyFigure.png」が設定されている。画像ファイルは2枚保存され、一枚は設定通りの名前、もう一枚はファイル名の最後が「_sgm.png」に変更されたものである。

保存先ファイル名を設定して図3の「保存」ボタンをクリックすると計算が始まるが、Stanスクリプトのコンパイルとビルドにかなりの時間が掛かる。ビルド後にMCMCサンプリングが始まり、サンプリングが終了すると平均の事後分布を描いたグラフが表示される(図4)。

図4

 

図4のウィンドウを閉じると、標準偏差の事後分布のグラフが表示される(図5)。

図5

 

図5のウィンドウを閉じるとプログラムの実行終了である。表示されたグラフは、ファイルに保存されている。

 

 

リスト1 1変量正規分布用Stanスクリプト(OneVarNormal.stan

 

data {

    int N;

    real X[N];

}

parameters {

    real mu;

    real<lower = 0.0> sgm;

}

model {

    mu ~ normal(0.0, 10000.0);

    sgm ~ exponential(0.0001);

    for (i in 1:N) {

        X[i] ~ normal(mu, sgm);

    }

}

 

 

リスト2 ベイズ分析Pythonスクリプト(MainOneVar.py

 

import tkinter as tk

import tkinter.filedialog as fd

import pandas as pd

import numpy as np

import pystan

import matplotlib.pyplot as plt

import seaborn as sb

 

def askinfilename(msg = None, types = [('', '*.*')]):

    """

            入力用ファイル名の設定

    """

    rt = tk.Tk()

    rt.withdraw()

    flnm = fd.askopenfilename(title = msg, filetypes = types)

    rt.destroy()

    return flnm

 

def askoutfilename(msg = None, types = [('', '*.*')], deftype = '.txt'):

    """

            出力用ファイル名の設定

    """

    rt = tk.Tk()

    rt.withdraw()

    flnm = fd.asksaveasfilename(title = msg, filetypes = types,

                                           defaultextension = deftype)

    rt.destroy()

    return flnm

 

 

finnm = askinfilename('入力データExcelファイル', types = [('Excel(*.xlsx)', '*.xlsx')])

foutnm = askoutfilename(msg = '画像保存ファイル', types = [('画像 (*.png)', '*.png')],

                       deftype = '.png')

 

xlsx_fl = pd.ExcelFile(finnm)

data_xlsx = pd.read_excel(xlsx_fl,  header = None)

data = data_xlsx.values

 

VarName = data[0,1]

X = data[1:, 1]

X = list(X)

print(X)

 

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

 

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

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

print(fit)

 

mu_025, mu_50, mu_975 = np.percentile(fit['mu'], [2.5, 50, 97.5])

 

#       日本語表記の設定

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

 

plt.title('平均値の事後分布' + \

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

              mu_50, mu_025, mu_975), fontsize = 20)

sb.kdeplot(fit['mu'])

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

plt.yticks([])

plt.tight_layout()

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

plt.show()

 

sgm_025, sgm_50, sgm_975 = np.percentile(fit['sgm'], [2.5, 50, 97.5])

plt.title('標準偏差の事後分布' + \

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

              sgm_50, sgm_025, sgm_975), fontsize = 20)

sb.kdeplot(fit['sgm'])

plt.xlabel(VarName + 'の標準偏差', fontsize = 14)

plt.yticks([])

plt.tight_layout()

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

plt.show()

 

 

Up