Edited at

# 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')
```

## (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')

np.random.seed(1234)
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)

```

## (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()

np.random.seed(1234)
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()
```