StanとRでベイズ統計モデリング(アヒル本)をPythonにしてみる - 7.2 対数をとるか否か

[実行環境](https://qiita.com/shngt/items/ce6f2c3ef348cd80f432)

### インポート

```python
import numpy as np
from scipy import stats
import pandas as pd
import pystan
import matplotlib.pyplot as plt
from matplotlib.figure import figaspect
from matplotlib.ticker import ScalarFormatter
import seaborn as sns
%matplotlib inline
```

### データ読み込み

```python
rental = pd.read_csv('./data/data-rental.txt')
```

# 7.2 対数をとるか否か

```python
_, (ax1, ax2) = plt.subplots(1, 2, figsize=figaspect(3/8))
rental.plot.scatter('Area', 'Y', ax=ax1)
plt.setp(ax1, xticks=np.arange(20, 120, 20), yticks=np.arange(500, 2000, 500))
rental.plot.scatter('Area', 'Y', ax=ax2)
xticks = (1,2,5,10,20,50,100)
yticks = (200,500,1000,2000)
ax2.get_xaxis().set_major_formatter(ScalarFormatter())
plt.setp(ax2, xscale='log', yscale='log', yticks=(10,20,50,100,200,500,1000,2000))
plt.setp(ax2, xticks=xticks, xticklabels=xticks, yticks=yticks, yticklabels=yticks, xlim=(9, 120))
plt.show()
```
![fig7-1.png](https://qiita-image-store.s3.amazonaws.com/0/103085/c6bd80b4-f3b7-0df1-9485-a1aefa2b9264.png)

```python
N_new = 50
Area_new = np.linspace(10, 120, N_new)
data = {col: rental[col] for col in rental.columns}
data['N'] = rental.index.size
data['N_new'] = N_new
data['Area_new'] = Area_new
fit1 = pystan.stan('./stan/model7-1.stan', data=data, seed=1234)
```

```python
for key in ['Area', 'Y', 'Area_new']:
    data[key] = np.log10(data[key])
fit2 = pystan.stan('./stan/model7-2.stan', data=data, seed=1234)
```

```python
ms1 = fit1.extract()
ms2 = fit2.extract()
for key in ['y_pred', 'y_new']:
    ms2[key] = np.power(10, ms2[key])
_, axes = plt.subplots(1, 2, figsize=figaspect(3/8), sharex=True, sharey=True)
for ms, ax in zip([ms1, ms2], axes):
    d_est = np.percentile(ms['y_new'], (10, 25, 50, 75, 90), axis=0)
    ax.fill_between(Area_new, d_est[0], d_est[-1], color='k', alpha=0.2)
    ax.fill_between(Area_new, d_est[1], d_est[-2], color='k', alpha=0.4)
    ax.plot(Area_new, d_est[2], color='k')
    rental.plot.scatter('Area', 'Y', ax=ax, color='k')
plt.show()
```
![fig7-2.png](https://qiita-image-store.s3.amazonaws.com/0/103085/75a1689f-2c3e-07ad-1e7c-ed505a906222.png)

```python
_, axes = plt.subplots(1, 2, figsize=figaspect(1/2), sharex=True, sharey=True)
for ms, ax in zip([ms1, ms2], axes):
    d_est = np.percentile(ms['y_pred'], (10, 50, 90), axis=0)
    ax.errorbar(rental['Y'], d_est[1], yerr=[d_est[1]-d_est[0], d_est[-1]-d_est[1]], fmt='o', color='k', elinewidth=1)
    xmin, xmax = ax.get_xlim()
    ymin, ymax = ax.get_ylim()
    lim = (min(xmin, ymin), max(xmax, ymax))
    ax.plot(lim, lim, c='k', linestyle='dashed')
    plt.setp(ax, xlim=lim, ylim=lim, xlabel='Obserbed', ylabel='Predicted')
plt.show()
```
![fig7-3.png](https://qiita-image-store.s3.amazonaws.com/0/103085/d2490e8b-1316-37f1-765b-401a4a730690.png)

```python
N_mcmc1 = len(ms1['lp__'])
N_mcmc2 = len(ms2['lp__'])
d_noise1 = pd.DataFrame(rental['Y'].values - ms1['mu'], columns=['noise{}'.format(i) for i in range(rental.index.size)])
d_noise2 = pd.DataFrame(np.log10(rental['Y'].values) - ms2['mu'], columns=['noise{}'.format(i) for i in range(rental.index.size)])
d_est1 = d_noise1.join(pd.Series(np.arange(N_mcmc1), name='mcmc'))
d_est2 = d_noise2.join(pd.Series(np.arange(N_mcmc2), name='mcmc'))

def get_map(col):
    kernel = stats.gaussian_kde(col)
    dens_x = np.linspace(col.min(), col.max(), kernel.n)
    dens_y = kernel.pdf(dens_x)
    mode_i = np.argmax(dens_y)
    mode_x = dens_x[mode_i]
    mode_y = dens_y[mode_i]
    return pd.Series([mode_x, mode_y], index=['X', 'Y'])

d_mode1 = d_noise1.apply(get_map).T
d_mode2 = d_noise2.apply(get_map).T

s_MAP1 = get_map(ms1['s_Y'])['X']
s_MAP2 = get_map(ms2['s_Y'])['X']

_, axes = plt.subplots(1, 2, figsize=figaspect(3/8))
for d_mode, s_MAP, ax in zip([d_mode1, d_mode2], [s_MAP1, s_MAP2], axes):
    sns.distplot(d_mode['X'], bins=16, color='k', hist_kws={'facecolor': 'w', 'edgecolor': 'k'}, kde_kws={'shade': True}, ax=ax)
    xdata = ax.lines[0].get_xdata()
    ax.plot(xdata, stats.norm.pdf(xdata, scale=s_MAP), color='k', linestyle='dashed')
    plt.setp(ax, xlabel='value')

plt.show()
```
![fig7-4.png](https://qiita-image-store.s3.amazonaws.com/0/103085/f4ccc0ef-2dfd-0062-9983-b716fd714df3.png)