StanとRでベイズ統計モデリング(アヒル本)をPythonにしてみる - 10.1 パラメータの識別可能性

15092 ワード
PyStan Stan Python Python テキストリンク

実行環境

インポート

import numpy as np
from scipy import stats
import pandas as pd
import statsmodels.api as sm
import pystan
import matplotlib.pyplot as plt
from matplotlib.figure import figaspect
%matplotlib inline

データ読み込み

category = pd.read_csv('./data/data-category.txt')
usagitokame = pd.read_csv('./data/data-usagitokame.txt')

10.1 パラメータの識別可能性

10.1.3 ラベルスイッチング

np.random.seed(123)

N = 200
a = 0.3
d1 = np.random.normal(loc=0, scale=2, size=N)
d2 = np.random.normal(loc=-4, scale=1, size=N)
d3 = np.random.uniform(size=N) <= a
Y = d1*d3 + d2*(1-d3)
X = np.linspace(Y.min(), Y.max(), 100)

plt.figure(figsize=figaspect(3/4))
ax = plt.axes()
ax.hist(Y, facecolor='w', edgecolor='k', density=True)
ax.plot(X, a*stats.norm.pdf(X, loc=0, scale=2), c='k', linestyle='dashed')
ax.plot(X, (1-a)*stats.norm.pdf(X, loc=-4, scale=1), c='k', linestyle='dashed')
ax.plot(X, a*stats.norm.pdf(X, loc=0, scale=2) + (1-a)*stats.norm.pdf(X, loc=-4, scale=1), c='k', alpha=0.4, linestyle='solid', linewidth=3)
plt.setp(ax, xlabel='y', ylabel='density')
plt.show()

10.1.4 多項ロジスティック回帰

X = sm.add_constant(category.iloc[:, :-1])
X['Age'] /= 100
X['Income'] /= 1000
data = dict(
    N=category.index.size,
    D=X.columns.size,
    K=category['Y'].max(),
    X=X,
    Y=category['Y']
)
fit = pystan.stan('./stan/model10-2.stan', data=data, seed=1234)

10.1.5 ウサギとカメ

np.random.seed(123)

G = 30
mu_pf = (0, 1.5)
pf = np.random.normal(loc=mu_pf, scale=1, size=(G, 2))
d = pd.DataFrame(np.argsort(pf, axis=1)+1, columns=('Loser', 'Winner'))
tbl = pd.crosstab([], d['Winner'])
tbl

data = dict(
    N=usagitokame.columns.size,
    G=usagitokame.index.size,
    LW=usagitokame
)
fit = pystan.stan('./stan/model10-3.stan', data=data, seed=1234)