「実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門」第１部 まとめ

* [実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門](https://www.amazon.co.jp/exec/obidos/ASIN/4065165369/hatena-blog-22/)を読んだので、復習も兼ねて第１部のメモ
* [RとStanではじめる　ベイズ統計モデリングによるデータ分析入門：サポートページ](https://logics-of-blue.com/r-stan-bayesian-model-intro-book-support/)
* 読んだときに作ったR、Stanのファイルは[ココ](https://github.com/kzy611/r_stan_intro)

## 何が書いてある本か

* 第１部：ベイズ統計モデリングによるデータ分析の方法。（**このページ**）
* 第２部～：RとStanを用いて実際にサンプルデータの可視化、分析、推定結果の図示を行う方法（**後日、別でまとめたい**）

## ベイズ統計モデリング

統計モデル・ベイズ推論によってデータを分析する手法。また、この本ではモデルや事後確率分布の評価にMCMC（マルコフ連鎖モンテカルロ法、Markov chain Monte Carlo methods）を用いている

## 統計モデル

* 統計モデル：データに適合するように構築された確率モデル
* 確率モデル：確率的な表現が用いられた数理モデル
* 数理モデル：数式を用いて表現されたモデル
* モデル：観測したデータを生み出す確率的な過程を簡潔に記述したもの
* モデリング：モデルを作る行為。

### 確率モデルで使われる確率分布

他にもいろいろあるが、代表的なもの

* 離散型：離散一様分布、ベルヌーイ分布、二項分布、ポアソン分布
* 連続型：連続一様分布、正規分布、対数正規分布、t分布

### 統計モデルの例

$f(Y)$は、$Y$が出る確率

* 二項分布：コインの表が出る確率を$p$として、10回投げて$Y$回表が出る

```math
Y \sim Binom(10, p) \\
f(Y|N, p) = _{10}C_Y p^Y (1-p)^{N-Y}
```

* 正規分布：平均$\mu$、標準偏差$\sigma$に従う商品の売上高$Y$

```math
Y \sim Normal(\mu, \sigma^2) \\
f(Y) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp(- \frac {(y-\mu)^2}{2 \sigma^2})
```

* 正規分布：平均$\mu$が気温$x$によって変わり、標準偏差$\sigma$に従う商品の売上高$Y$

```math
\mu = \beta_0 + \beta_1 x \\
Y \sim Normal(\mu, \sigma^2) \\
f(Y) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp(- \frac {(y-\mu)^2}{2 \sigma^2})
```

## ベイズ推論

ベイズの定理を用いたベイズ更新によって、興味の対象となる条件付き確率を得ること

### ベイズの定理

$A$が発生する確率を$P(A)$、$X$が発生したという条件で$A$が発生する確率を$P(A|X)$とすると

```math
P(X|A) = \frac {P(A|X) P(X)} {P(A)}
```

$P(A) \neq 0$の条件で以下の式を変形した、と考えるとわかりやすい。（両辺とも$A$と$X$が同時に発生する確率）

```math
P(X|A) P(A) = P(A|X) P(X)
```

### ベイズ更新

ベイズの定理の$A$を$D$にして少し書き方を変えると

```math
P(X|D) = \frac {P(D|X)} {P(D)} P(X)
```

となる。このとき、$X$を確率分布のパラメータ（統計モデルの$p$とか$\mu$とか$\sigma$）、$D$を実際に得られたデータとすると、この式の意味は

「事前確率（何もデータがないときのパラメータ$X$の確率）に、データ$D$の情報を使うと、事後確率（データ$D$が得られたという条件でのパラメータ$X$の確率）に更新される」

という意味になる。

### 例題（モンティ・ホール問題）

* A, B, Cのドアのうち1つが正解の場合、事前確率は$P(A) = P(B) = P(C) = \frac{1}{3}$
* 回答者がドアAを選び、出題者がドアCを開けて「ドアCは正解ではない」と示す（この情報をDとする）
* このとき、$P(D), P(D|A), P(D|B), P(D|C)$は「回答者がドアAを選ぶ確率×出題者がドアCを開ける確率」なので以下のようになる。（出題者は正解を知ってるので、回答者が選んだドアと正解のドアを避けて開く）

```math
P(D) = \frac{1}{3} \frac{1}{2} = \frac{1}{6} \\
P(D|A) = \frac{1}{3} \frac{1}{2} = \frac{1}{6} \\
P(D|B) = \frac{1}{3} 1 = \frac{1}{3} \\
P(D|C) = \frac{1}{3} 0 = 0 \\
```

* これらを使ってベイズ更新すると、$P(A|D), P(B|D), P(C|D)$は以下になる

```math
P(A|D) = \frac {P(D|A)} {P(D)} P(A) = \frac {\frac{1}{6}} {\frac{1}{6}} \frac{1}{3} = \frac{1}{3} \\
P(B|D) = \frac {P(D|B)} {P(D)} P(B) = \frac {\frac{1}{3}} {\frac{1}{6}} \frac{1}{3} = \frac{2}{3} \\
P(C|D) = \frac {P(D|C)} {P(D)} P(C) = \frac {0} {\frac{1}{6}} \frac{1}{3} = 0 \\
```

### ベイズ統計の解説動画

ベイズ統計については、この動画がわかり易かった
[AIcia Solid Project - ベイズ統計](https://www.youtube.com/playlist?list=PLhDAH9aTfnxIU4Hd1G1UdIVzHpgKfyEnw)

## MCMC（マルコフ連鎖モンテカルロ）法

* MCMC（Markov chain Monte Carlo methods）法：マルコフ連鎖を利用して乱数を生成する手法。
* 上の例題では簡単に事後確率を求められたが、実際には求めるのに複雑な計算が必要な場合が多い。（特に、パラメータやデータが連続型だと複雑な積分計算が必要）
* そのため、頑張って全部は計算しないで、MCMCで事後確率分布に従う乱数を生成して、その乱数の平均値などを使う

### モンテカルロ法

* モンテカルロ法（Monte Carlo method）：シミュレーションや数値計算を乱数を用いて行う手法の総称
* 今回は、積分計算を乱数を使って代替する（モンテカルロ積分）

#### 例

求めたいパラメータ$\theta$があり、データ$D$が得られたときの事後確率分布を$f(\theta|D)$とすると、$\theta$の期待値は以下になる。（「$\theta$が発生する確率と$\theta$をかけあわせたもの」を全パターン足したもの）

```math
\int_{-\infty}^{\infty}f(\theta|D)d\theta
```

このとき$f(\theta|D)$が複雑すぎると積分計算できない。

ここで、かわりに事後分布に従う乱数$\hat{\theta}$を1000個（$\hat{\theta}_1, \hat{\theta}_2, ..., \hat{\theta} _{1000}$）発生できたとすると、期待値は以下のようになる。

```math
\frac{\sum_{i=1}^{1000} \hat{\theta_i}}{1000}
```

厳密な答えとは多少ズレるが、乱数の数を十分に大きくすれば実用的には問題ない近い値になる。

### マルコフ連鎖

* あるステップでの状態が、その直前の状態のみに依存するモデル（２ステップ以上前には依存しない）
* 例：今日の天気は、昨日の天気によって決まる。（一昨日の天気は関係ない）

#### 例題

以下の場合、全体でどのくらい晴れ・曇り・雨になるのかを求める

* 今日が晴れだったら、次の日は晴れ70%、曇り20%、雨10%
* 今日が曇りだったら、次の日は晴れ30%、曇り50%、雨20%
* 今日が雨だったら、次の日は晴れ20%、曇り40%、雨40%

[パターン1]

* 150個の場所でスタート時点ですべて雨として、200回更新する
* 晴れ、曇り、雨のどれも、ある程度の確率で安定している

![1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/563901/2780adcf-bd8f-c2e0-dc99-f1f0338e22ad.png)

[パターン2]

* 150個の場所でスタート時点ですべて晴れとして、10回更新する
* コレだと、収束してるのかどうかわからない = 回数が少なすぎると収束しない

![2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/563901/90557a01-66af-86d7-b5a5-51f8cdc7bc7f.png)

```python
# ライブラリ読み込み&初期設定
import numpy as np
import pandas as pd
import random

import matplotlib.pyplot as plt
%matplotlib inline

# 初期の配分
scr_cnt = {
    's': 150,
    'c': 0,
    'r': 0,
}

# 繰り返し回数
iter_num = 10
# 乱数設定
seed = 0
```

```python
# 天気更新用の関数
# 晴れs:0、くもりc:1、雨r:2
def update_w(w):
    rand = random.random()
    
    # 移行確率
    ratio = {
        0: {0: 0.7, 1: 0.2, 2: 0.1}, # 晴れ  - >晴れ、くもり、雨
        1: {0: 0.3, 1: 0.5, 2: 0.2}, # くもり-> 晴れ、くもり、雨
        2: {0: 0.2, 1: 0.4, 2: 0.4}, # 雨-    > 晴れ、くもり、雨
    }
    
    r_cum = 0
    for i, r in ratio[w].items():
        r_cum += r
        if rand < r_cum:
            return i
```

```python
# 乱数を生成
# 初期状態
w = np.concatenate((np.repeat(0, scr_cnt['s']),
                    np.repeat(1, scr_cnt['c']),
                    np.repeat(2, scr_cnt['r'])),
                   axis=0)

# 結果格納用
x = []
s_ratios = []
c_ratios = []
r_ratios = []

# 初期値格納
s_ratio = sum(x==0 for x in w) / len(w)
c_ratio = sum(x==1 for x in w) / len(w)
r_ratio = sum(x==2 for x in w) / len(w)
x.append(0)
s_ratios.append(s_ratio)
c_ratios.append(c_ratio)
r_ratios.append(r_ratio)

random.seed(seed)
for i in range(iter_num):
    w = list(map(update_w, w)) # 天気更新
    
    s_ratio = sum(x==0 for x in w) / len(w)
    c_ratio = sum(x==1 for x in w) / len(w)
    r_ratio = sum(x==2 for x in w) / len(w)
    x.append(i+1)
    s_ratios.append(s_ratio)
    c_ratios.append(c_ratio)
    r_ratios.append(r_ratio)
```

```python
# プロット
fig = plt.figure()
ax = fig.add_subplot(111)

ax.plot(x, s_ratios, label='sunny')
ax.plot(x, c_ratios, label='cloudy')
ax.plot(x, r_ratios, label='rainy')
ax.set_xlim([0, iter_num])
ax.set_ylim([0, 1])
ax.legend(loc='best')

plt.show()
```

```python
# 乱数の統計量を算出
df_ratios = pd.DataFrame({
    'sunny': s_ratios,
    'cloudy': c_ratios,
    'rainy': r_ratios,
})

# 初期値の影響を受けやすい値(今回は前半の半分)を除外
df_ratios.loc[len(df_ratios)//2:, :].describe()
```

### 事後分布に従う乱数の生成

* この本では「メトロポリス・ヘイスティングス法（MH法）」と「ハミルトニアンモンテカルロ法（HMC法）」について書かれている
* ここにはMH法のみ書く（HMC法は、都度ググる）

#### メトロポリス・ヘイスティングス法（MH法）

1. データ$D$が得られたとする
2. データが生成されるモデル（確率分布）を考え、適当にパラメータの初期値$\hat{\theta}_1$を決める
3. $\hat{\theta_1}$に適当な乱数（正規分布で生成する）を足して、$\hat{\theta_2}^{仮}$とする
4. データ$D$が得られた状態での事後確率$f(\hat{\theta_1}|D), f(\hat{\theta_2}^{仮}|D)$を比べ、決められたルールに従い$\hat{\theta_2}$を決定する
5. 同様に、$i=3, 4, \cdots$についても$\hat{\theta_i}$を求めていく
6. 一定数$\hat{\theta}$が得られたら、初期値の影響を受ける最初の方の$\hat{\theta}$を除外し、残った$\hat{\theta}$の平均値などをとってパラメータの点推定値とする

**$\hat{\theta_i}$の決定方法**

ベイズの定理より

```math
f(\theta|D) = \frac {f(D|\theta) f(\theta)} {f(D)}
```

なので、$\hat{\theta_2}^{仮}, \hat{\theta_1}$の事後確率を比較した$rate$は、$f(D|\theta) f(\theta) = Kernel(\theta)$とおくと


```math
\begin{align}
rate &= \frac {f(\hat{\theta_2}^{仮}|D)}{f(\hat{\theta_1}|D)} \\
&= \frac {f(D|\hat{\theta_2}^{仮}) f(\hat{\theta_2}^{仮})} {f(D)} \frac {f(D)} {f(D|\hat{\theta_1}) f(\hat{\theta_1})} \\
&= \frac {f(D|\hat{\theta_2}^{仮}) f(\hat{\theta_2}^{仮})} {f(D|\hat{\theta_1}) f(\hat{\theta_1})} \\
&= \frac {Kernel(\hat{\theta_2}^{仮})} {Kernel(\hat{\theta_1})} \\
\end{align}
```

* rateが1を超えると、$\hat{\theta_2}^{仮}$のほうが事後確率が高い
* 計算するのが難しい$f(D)$が消えており、計算可能となっている

$rate$を求めた後は、以下のように$\hat{\theta_i}$を更新する

* $rate > 1$の場合：$\hat{\theta_i} = \hat{\theta_i}^{仮}$とする
* $rate < 1$の場合：確率$rate$で$\hat{\theta_i} = \hat{\theta_i}^{仮}$、確率$(1-rate)$で$\hat{\theta_i} = \hat{\theta_{i-1}}$とする

「変更したほうが事後確率が上がる場合、確定で変更」「変更したほうが事後確率が下がる場合も、変更することもある」

#### 得られた乱数の取り扱いの注意

* 各種設定：chains、iter、warmup、thin
* 収束の評価
* 代表値のとり方：事後中央値、事後期待値（平均）、事後確率最大値（最頻値）

## 他参考リンク

* [Stan超初心者入門](https://www.slideshare.net/simizu706/stan-62042940)
* [確率分布（Wikipedia）](https://ja.wikipedia.org/wiki/%E7%A2%BA%E7%8E%87%E5%88%86%E5%B8%83)
* [いろいろな確率分布の平均，分散，特性関数などまとめ](https://mathtrain.jp/distributionmatome)
* [ベイズ推定（Wikipedia）](https://ja.wikipedia.org/wiki/%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E5%AE%9A)
* [理由不十分の原則](https://clover.fcg.world/2016/07/30/6249/)
* [モンティ・ホール問題](https://ja.wikipedia.org/wiki/%E3%83%A2%E3%83%B3%E3%83%86%E3%82%A3%E3%83%BB%E3%83%9B%E3%83%BC%E3%83%AB%E5%95%8F%E9%A1%8C)
* [モンテカルロ法（Wikipedia）](https://ja.wikipedia.org/wiki/%E3%83%A2%E3%83%B3%E3%83%86%E3%82%AB%E3%83%AB%E3%83%AD%E6%B3%95)
* マルコフ連鎖
  * [マルコフ連鎖の基本とコルモゴロフ方程式](https://mathtrain.jp/markovchain)
  * [Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%AB%E3%82%B3%E3%83%95%E9%80%A3%E9%8E%96)
* [メトロポリス・ヘイスティングス法（Wikipedia）](https://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%88%E3%83%AD%E3%83%9D%E3%83%AA%E3%82%B9%E3%83%BB%E3%83%98%E3%82%A4%E3%82%B9%E3%83%86%E3%82%A3%E3%83%B3%E3%82%B0%E3%82%B9%E6%B3%95)
* ハミルトニアンモンテカルロ法
  * [ハミルトニアンモンテカルロ法（Qiita）](https://qiita.com/hiro5585/items/f50df506b208d37cdd7d)
  * [ハミルトニアンモンテカルロ法（Wikipedia）](https://ja.wikipedia.org/wiki/%E3%83%8F%E3%83%9F%E3%83%AB%E3%83%88%E3%83%8B%E3%82%A2%E3%83%B3%E3%83%BB%E3%83%A2%E3%83%B3%E3%83%86%E3%82%AB%E3%83%AB%E3%83%AD%E6%B3%95)