「実践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)

## 何が書いてある本か

* 第１部：ベイズ統計モデリングによるデータ分析の方法。（[コッチのページ](https://qiita.com/kazuya_minakuchi/items/bfe90e62a11eda162956)）
* 第２部：RとStanの基本的な使い方
* 第３部～：サンプルデータを用いた、モデルの作り方（**このページ**）

## 一般化線形モデル（GLM: Generalized Linear Models）

* 確率分布、線形予測子、リンク関数の3つの要素からなるモデル

```math
\begin{align}
g(\mu_i) &= \beta_0 + \beta_1 x_{i} \\
y_i &\sim Poiss(\lambda_i) \\
\end{align}
```

とすると、

* 確率分布：$Poiss$
* 線形予測子：$\beta_0 + \beta_1 x_{i}$
* リンク関数：$g(\mu_i)$

### 例

* 平均$\mu$・分散$\sigma^2$の正規分布に従う$y$

```math
\begin{align}
y &\sim Normal(\mu, \sigma^2) \\
\end{align}
```

* $\mu_i$が$x_i$によって変化する、平均$\mu_i$・分散$\sigma^2$の正規分布に従う$y_i$

```math
\begin{align}
\mu_i &= \beta_0 + \beta_1 x_i \\
y_i &\sim Normal(\mu_i, \sigma^2) \\
\end{align}
```

* $\lambda_i$が$x_{i1}, x_{i2}, x_{i3}$によって変化する、ポアソン分布に従う$y_i$

```math
\begin{align}
log(\lambda_i) &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} \\
y_i &\sim Poiss(\lambda_i) \\
\end{align}
```

* $p_i$が$x_{i1}, x_{i2}$によって変化する、試行回数nの二項分布に従う$y_i$

```math
\begin{align}
logit(p_i) &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} \\
y_i &\sim Binom(n, p_i) \\
\end{align}
```

または

```math
\begin{align}
p_i &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} \\
y_i &\sim Binom(n, logistic(p_i)) \\
\end{align}
```

ここで

```math
\begin{align}
logit(p) &= log( \frac{p}{1-p} ) \\
logistic(x) &= \frac{1}{1+exp(-x)} \\
\end{align}
```

* $\mu_i$が$x_{i1}, x_{i2}$によって変化する、平均$\mu_i$・分散$\sigma^2$の正規分布に従う$y_i$。ただし、説明変数同時に交互作用も考慮する

```math
\begin{align}
\mu_i &= \beta_0 + \beta_{11} x_{i1} + \beta_{22} x_{i2} + \beta_{12} x_{i1} x_{i2} \\
y_i &\sim Normal(\mu_i, \sigma^2) \\
\end{align}
```

## 一般化線形混合モデル（GLMM: Generalized Linear Mixed Models）

* 一般化線形モデル（GLM）にランダム効果を加えて混合モデルとしたもの

### 例

* $\lambda_i$が$x_{i1}, x_{i2}$によって変化する、ポアソン分布に従う$y_i$。データごとに切片へのランダムな影響$r_i$を加える

```math
\begin{align}
r_i &\sim Normal(0, \sigma_r^2) \\
log(\lambda_i) &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_i \\
y_i &\sim Poiss(\lambda_i) \\
\end{align}
```

* $\lambda_i$が$x_{i1}, x_{i2}$によって変化する、ポアソン分布に従う$y_i$。データのグループごとに切片へのランダムな影響$r_k$を加える

```math
\begin{align}
r_k &\sim Normal(0, \sigma_r^2) \\
log(\lambda_i) &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_k \\
y_i &\sim Poiss(\lambda_i) \\
\end{align}
```

* $\lambda_i$が$x_{i1}$によって変化する、ポアソン分布に従う$y_i$。データのグループごとに切片、係数へのランダムな影響$r_k, \tau_k$を加える

```math
\begin{align}
r_k &\sim Normal(0, \sigma_r^2) \\
\tau_k &\sim Normal(0, \sigma_{\tau}^2) \\
log(\lambda_i) &= \beta_0 + (\beta_1 + \tau_k) x_{i1} + r_k \\
y_i &\sim Poiss(\lambda_i) \\
\end{align}
```

## 状態空間モデル

時系列データを、目に見えない状態を表す状態方程式と、そこに観測誤差などを加えて実際に計測されたデータを表す観測方程式に分けて考える

### 動的線形モデル（DLM: Dynamic Linear Models）

#### ローカルレベルモデル

状態$\mu$が「1つ前の時点+正規ホワイトノイズ」のランダムウォークし、それに正規分布にしたがうノイズが加わって観測値$y_t$が得られるモデル

```math
\begin{align}
\mu_t &= \mu_{t-1} + w_t, \quad w_t \sim Normal(0, \sigma_w^2) \\
y_t &= \mu_t + v_t, \quad v_t \sim Normal(0, \sigma_v^2) \\
\end{align}
```

書き換えると

```math
\begin{align}
\mu_t &\sim Normal(\mu_{t-1}, \sigma_w^2) \\
y_t &\sim Normal(\mu_t, \sigma_v^2) \\
\end{align}
```

#### 時変係数モデル

「ランダムウォークする切片」を持つローカルレベルモデルに、「ランダムウォークする回帰係数」を持つ説明変数を加えたモデル

```math
\begin{align}
\mu_t &\sim Normal(\mu_{t-1}, \sigma_w^2) \\
\beta_t &\sim Normal(\beta_{t-1}, \sigma_{\tau}^2) \\
\alpha_t &= \mu_t + \beta_t {ex}_t \\
y_t &\sim Normal(\alpha_t, \sigma_v^2) \\
\end{align}
```

#### 平滑化トレンドモデル

状態$\mu$の変化が、「1つ前の状態$\mu$の変化+正規ホワイトノイズ」に従うモデル

```math
\begin{align}
\mu_t - \mu_{t-1} &= \mu_{t-1} - \mu_{t-2} + \zeta_t, \quad \zeta_t \sim Normal(0, \sigma_{\zeta}^2) \\
\end{align}
```

を変形し、状態方程式と観測方程式は以下の形になる

```math
\begin{align}
\mu_t &= 2 \mu_{t-1} - \mu_{t-2} + \zeta_t, \quad \zeta_t \sim Normal(0, \sigma_{\zeta}^2) \\
y_t &\sim Normal(\mu_t, \sigma_v^2) \\
\end{align}
```

書き換えると

```math
\begin{align}
\mu_t &\sim Normal(2 \mu_{t-1} - \mu_{t-2}, \sigma_{\zeta}^2) \\
y_t &\sim Normal(\mu_t, \sigma_v^2) \\
\end{align}
```

#### ローカル線形トレンドモデル

平滑化トレンドモデルに、過程誤差を加えたもの。
平滑化トレンドモデルの式で

```math
\begin{align}
\delta_t &= \mu_t - \mu_{t-1}
\end{align}
```

と置くことで、平滑化トレンドモデルは

```math
\begin{align}
\delta_t &= \delta_{t-1} + \zeta_t, \quad \zeta_t \sim Normal(0, \sigma_{\zeta}^2) \\
\mu_t &= \mu_{t-1} + (\mu_{t-1} - \mu_{t-2} + \zeta_t) \\
 &= \mu_{t-1} + \delta_{t-1} \\
y_t &\sim Normal(\mu_t, \sigma_v^2) \\
\end{align}
```

となり、これに過程誤差$w_t$を加えたローカル線形トレンドモデルは以下の形になる

```math
\begin{align}
\delta_t &= \delta_{t-1} + \zeta_t, \quad \zeta_t \sim Normal(0, \sigma_{\zeta}^2) \\
\mu_t &= \mu_{t-1} + \delta_{t-1} + w_t, \quad w_t \sim Normal(0, \sigma_w^2) \\
y_t &\sim Normal(\mu_t, \sigma_v^2) \\
\end{align}
```

書き換えると

```math
\begin{align}
\delta_t &\sim Normal(\delta_{t-1}, \sigma_{\zeta}^2) \\
\mu_t &\sim Normal(\mu_{t-1} + \delta_{t-1}, \sigma_w^2) \\
y_t &\sim Normal(\mu_t, \sigma_v^2) \\
\end{align}
```

#### 周期性

周期性の成分は、「周期成分の合計値が正規ホワイトノイズ$s_t$に従う」というルールを適用する。

```math
\begin{align}
\gamma_1 + \gamma_2 + \cdots + \gamma_k &= \sum_{i=1}^{k} \gamma_i = s_t, \quad s_t \sim Normal(0, \sigma_s^2) \\
\end{align}
```

四半期ごとのデータなら$k=4$、日データなら$k=4$、月データなら$k=12$にする。これを変形すると以下の形になる。

```math
\begin{align}
\gamma_k + \sum_{i=1}^{k-1} \gamma_i &\sim Normal(0, \sigma_s^2) \\
\gamma_k &\sim Normal(- \sum_{i=1}^{k-1} \gamma_i, \sigma_s^2) \\
\end{align}
```

よってt時点での周期成分$gamma_t$は

```math
\begin{align}
\gamma_t &\sim Normal(- \sum_{i=t-(k-1)}^{t-1} \gamma_i, \sigma_s^2) \\
\end{align}
```


この成分を入れて「トレンド成分+周期成分+ホワイトノイズ」で表現される**基本構造時系列モデル（Basic Structual Time Series Modes）**の状態方程式と観測方程式は以下になる。

* 1次トレンド（ローカルレベルモデル）+周期成分

```math
\begin{align}
\mu_t &\sim Normal(\mu_{t-1}, \sigma_w^2) \\
\gamma_t &\sim Normal(- \sum_{i=t-(k-1)}^{t-1} \gamma_i, \sigma_s^2) \\
\alpha_t &= \mu_t + \gamma_t \\
y_t &\sim Normal(\alpha_t, \sigma_v^2) \\
\end{align}
```

* 2次トレンド（平滑化トレンドモデル）+周期成分

```math
\begin{align}
\mu_t &\sim Normal(2 \mu_{t-1} + \mu_{t-2}, \sigma_w^2) \\
\gamma_t &\sim Normal(- \sum_{i=t-(k-1)}^{t-1} \gamma_i, \sigma_s^2) \\
\alpha_t &= \mu_t + \gamma_t \\
y_t &\sim Normal(\alpha_t, \sigma_v^2) \\
\end{align}
```

#### 自己回帰モデル

過去の自分の値を説明変数に持つモデル

* 1次の自己回帰モデル

```math
\begin{align}
\alpha_t &\sim Normal(\beta_0 + \beta_1 \alpha_{t-1}, \sigma_w^2) \\
y_t &\sim Normal(\alpha_t, \sigma_v^2) \\
\end{align}
```

### 動的一般化線形モデル（DGLM: Dynamic Generalized Linear Models）

動的線形モデル（DLM: Dynamic Linear Models）のリンク関数や確率分布を変えて、正規分布以外の使えるようにしたもの

#### 二項分布

```math
\begin{align}
\mu_t &\sim Normal(\mu_{t-1}, \sigma_w^2) \\
y_t &\sim Bernoulli(logistic(\mu_t)) \\
\end{align}
```

#### ポアソン分布

* 2次トレンド+説明変数

```math
\begin{align}
\mu_t &\sim Normal(2 \mu_{t-1} - \mu_{t-2}, \sigma_{\zeta}^2) \\
\lambda_t &= \mu_t + \beta {ex}_t \\
y_t &\sim Poisson(exp(\lambda_t)) \\
\end{align}
```

* 2次トレンド+説明変数+ランダム効果

```math
\begin{align}
r_t &\sim Normal(0, \sigma_r^2) \\
\mu_t &\sim Normal(2 \mu_{t-1} - \mu_{t-2}, \sigma_{\zeta}^2) \\
\lambda_t &= \mu_t + \beta {ex}_t + r_t \\
y_t &\sim Poisson(exp(\lambda_t)) \\
\end{align}
```

## 他参考リンク

* [一般化線形モデル（Wikipedia）](https://ja.wikipedia.org/wiki/%E4%B8%80%E8%88%AC%E5%8C%96%E7%B7%9A%E5%BD%A2%E3%83%A2%E3%83%87%E3%83%AB)
* [30分だけでは決してよくわからないとてもとても難しい一般化線形モデル with R](http://cogpsy.educ.kyoto-u.ac.jp/personal/Kusumi/datasem13/shrasuna1.pdf)
* 状態空間モデル
  * [状態空間モデル](https://logics-of-blue.com/%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AB%E3%81%A8%E3%81%AF/)
  * [なぜ状態空間モデルを使うのか](https://logics-of-blue.com/%E3%81%AA%E3%81%9C%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AB%E3%82%92%E4%BD%BF%E3%81%86%E3%81%AE%E3%81%8B/)