StanとRでベイズ統計モデリング(アヒル本)をPythonにしてみる - Chapter 12 練習問題

実行環境

インポート

import numpy as np
import pandas as pd
import pystan
import matplotlib.pyplot as plt
from matplotlib.figure import figaspect
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
%matplotlib inline

データ読み込み

ss2 = pd.read_csv('./data/data-ss2.txt')
mesh2d = pd.read_csv('./data/data-2Dmesh.txt', header=None)
mesh2d_design = pd.read_csv('./data/data-2Dmesh-design.txt', header=None)

(1)

np.random.seed(123)

def SimulateSS(N, T, s_mu, s_Y):
    d = pd.DataFrame()
    for n in range(N):
        mu = np.zeros(T)
        mu[0] = 10
        for t in range(1, T):
            mu[t] = np.random.normal(mu[t-1], s_mu)
        Y = np.random.normal(mu, s_Y)
        d = d.append(pd.DataFrame(dict(Trial=n+1, Time=np.arange(T)+1, Y=np.round(Y, 2))), ignore_index=True)
    return d

d1 = SimulateSS(N=5, T=21, s_mu=2, s_Y=0.1)
d2 = SimulateSS(N=5, T=21, s_mu=0.1, s_Y=2)

_, axes = plt.subplots(1, 2, figsize=figaspect(1/2))
for d, ax in zip([d1, d2], axes):
    sns.lineplot('Time', 'Y', hue='Trial', data=d, ci=None, ax=ax)
plt.show()

(2)

data = dict(
    T=ss2.index.size,
    T_pred=8,
    Y=ss2['Y']
)
stanmodel = pystan.StanModel('./stan/chap12ex2.stan')
fit = stanmodel.sampling(data=data, iter=4000, thin=5, seed=1234)

(3)

lowess用にRを使用するための準備。

import rpy2.robjects as ro
from rpy2.robjects import numpy2ri
numpy2ri.activate()

r_command = '''
I <- nrow(d)
J <- ncol(d)
rownames(d) <- 1:I
colnames(d) <- 1:J
d_melt <- reshape2::melt(d)
colnames(d_melt) <- c('i','j','Y')
loess_res <- loess(Y ~ i + j, data=d_melt, span=0.1)
smoothed <- matrix(loess_res$fitted, nrow=I, ncol=J)
'''

d_ori = mesh2d
TID = mesh2d_design
T = TID.values.max()

stanmodel = pystan.StanModel('./stan/model12-13.stan')

def estimate(s_add):
    np.random.seed(1234)
    d = np.random.normal(loc=d_ori, scale=s_add)
    ro.r.assign('d', d)
    ro.r(r_command)

    data = dict(
        I=d_ori.index.size,
        J=d_ori.columns.size,
        Y=d,
        T=T,
        TID=TID
    )
    fit = stanmodel.sampling(data=data, iter=5000, thin=5, seed=1234, init=lambda: dict(
        r=ro.r('smoothed'),
        s_r=1,
        s_Y=1,
        s_beta=1,
        beta=np.random.normal(0, 0.1, T)
    ), n_jobs=-1)

res = [estimate(s_add) for s_add in (0.1, 0.2, 0.3)]

(4)

r_command = '''
I <- nrow(d)
J <- ncol(d)
rownames(d) <- 1:I
colnames(d) <- 1:J
d_melt <- reshape2::melt(d)
colnames(d_melt) <- c('i','j','Y')
loess_res <- loess(Y ~ i + j, data=d_melt, span=0.1)
smoothed <- matrix(loess_res$fitted, nrow=I, ncol=J)
'''

d_ori = mesh2d
I=d_ori.index.size
J=d_ori.columns.size
TID = mesh2d_design
T = TID.values.max()
s_add = 0.3

np.random.seed(1234)
d = np.random.normal(loc=d_ori, scale=s_add)
ro.r.assign('d', d)
ro.r(r_command)

stanmodel = pystan.StanModel('./stan/chap12ex4.stan')

data = dict(
    I=I,
    J=J,
    Y=d,
    T=T,
    TID=TID,
    S_s_Y=0.1
)
fit = stanmodel.sampling(data=data, iter=5000, thin=5, seed=1234, init=lambda: dict(
    r=ro.r('smoothed'),
    s_r=1,
    s_Y=1,
    s_beta=1,
    beta=np.random.normal(0, 0.1, T)
))

%matplotlib notebook
%matplotlib notebook

ms = fit.extract()

r_median = np.median(ms['r'], axis=0)
ii, jj = np.mgrid[:I, :J]

ax = Axes3D(plt.gcf())
ax.plot_wireframe(ii, jj, r_median, color='k')
ax.plot_surface(ii, jj, r_median, color='k', alpha=0.2)
plt.setp(ax, xlabel='Plate Row', ylabel='Plate Column', zlabel='r')
plt.show()