---
tags: Survey
lang: ja-jp
---
# Gibbs Sampling & Change Point Detection

## Next Action
- 「周期性を考慮した」変化点検出の手法の調査

## Gibbs Sampling
> Griffiths & Steyvers 2004

### TL;DR
1. ベイズ推論の目標は**確率変数の集合について完全な事後確率分布を獲得する**こと
2. Gibbs Samplingは確率変数の事後分布を獲得するための**MCMC**な手法の一つ
3. Gibbs Samplingでは，**残りの変数を現在の値に固定した条件付き分布から掃き出す**ことで，各確率変数（のかたまり）の事後分布を獲得する
4. MCMCの理論は，Gibbs Samplingを通じて得られる**確率変数の事後分布は各確率変数の同時事後分布に収束する**ことを保証する

### MCMC Samplingの背景
MCMCの考え方は，ある確率分布に従う確率変数の期待値はエルゴード的な平均から推定することができるというもの．$\mathcal{P}$を事後分布とし，$s^{(i)}$をその確率変数とすると，$f$による像の期待値は次のような式で推定することができる．

$$
E[f(s)]_{\mathcal{P}} \approx \frac{1}{N}\sum_{i=1}^{N}f(s^{(i)})
$$

この式は事後分布が具体的に分かっていなくても，**$N$個のサンプルから導出される代表値を確率変数の期待値として採用可能である**ことを意味している．

### Gibbs Sampling Algorithm
「サンプリングの結果から期待値を推定する」というのがMCMCの基本的な考え方．
ここでの問題は，実際にどのようにして事後分布からサンプルを得るのかの方法についてである．Gibbs Samplingは事後分布（から生成される確率変数）を，残りの確率変数を現在の値に固定した条件付き分布から生成されたものから求める手法である．

#### Algorithm
- $x^{(0)}$を事前分布$q$からサンプリングし，それを初期値とする．
- $i=1,\ldots,M$について，以下を繰り返す．
- $x_{1}^{(i)} \sim p\left(X_{1}=x_{1} | X_{2}=x_{2}^{(i-1)}, X_{3}=x_{3}^{(i-1)}, \ldots, X_{D}=x_{D}^{(i-1)}\right)$
- $x_{2}^{(i)} \sim p\left(X_{2}=x_{2} | X_{1}=x_{1}^{(i)}, X_{3}=x_{3}^{(i-1)}, \ldots, X_{D}=x_{D}^{(i-1)}\right)$
- $x_{D}^{(i)} \sim p\left(X_{D}=x_{D} | X_{1}=x_{1}^{(i)}, X_{2}=x_{2}^{(i)}, \ldots, X_{D}=x_{D-1}^{(i)}\right)$

:::info
$x$は確率変数の集合（変数を束ねたベクトル）とし，右上の添字はサンプリングの回数，右下の添字は確率変数のラベル（何次元目か）を表現する．注意すべきは，$i$番目の次元の確率変数をサンプリングした結果を利用して$i+1$番目の次元の確率変数のサンプリングを行っていること．
:::

$M$の値は通常大きなものが設定される．イテレーションの初期にサンプリングされた結果は事後分布に従わないことが通常であるので利用されない．この利用しないイテレーションの数を`burn-in`と呼んでいる．

## Change Point Detection
時系列の観測データから，傾向が変わった時点を検出するための手法．

### カウントデータの変化点検出
#### 問題設定
$x_{1}, x_{2}, \ldots, x_{N}$で表現される時系列のカウントデータを考える．
このデータは$1,\ldots,n$と$n+1,\ldots,N$の２つの区間でカウントの平均値が異なる．
ここでは，$i$回目の観測値はポアソン分布に従う確率変数としてモデリングする．

ポアソン分布の確率関数は次の通り：

$$
\begin{aligned}
\operatorname{Poisson}(x ; \lambda)
&=e^{-\lambda} \frac{\lambda^{x}}{x !} \\
&=\exp (x \log \lambda-\lambda-\log (x !)) \end{aligned}
$$

$\lambda$はポアソン分布のパラメータであり，平均を意味する．この平均値$\lambda$をガンマ分布に従って生成される確率変数としてモデリングする．

ガンマ分布の確率（密度）関数は次の通り：

$$
\begin{aligned} \operatorname{Gamma}(\lambda ; a, b)
&=\frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) \\
&=\exp ((a-1) \log \lambda-b \lambda-\log \Gamma(a)+a \log b) \end{aligned}
$$

$n$回目の観測値が従うポアソン分布のパラメータを$\lambda_{1}$とし，それ以降は$\lambda_{2}$にパラメータが変化するという問題を扱っていることになる．

生成モデルの過程を整理すると次のようになる：

$$
\begin{aligned} n & \sim \text { Uniform }(1,2, \ldots, N) \\
\lambda_{i} & \sim \operatorname{Gamma}\left(\lambda_{i} ; a, b\right) \\
x_{i} & \sim\left\{\begin{array}{ll}{\operatorname{Poisson}\left(\lambda_{i} ; \lambda_{1}\right)} & {1 \leq i \leq n}
\\ {\text { Poisson }\left(x_{i} ; \lambda_{2}\right)} & {n<i \leq N}\end{array}\right.\end{aligned}
$$

潜在変数$n$，$\lambda_{1}$，$\lambda_{2}$の事後分布の推論問題は（同時分布からベイズの定理を利用することで）次のように整理できる：

$$
p\left(\lambda_{1}, \lambda_{2}, n | x_{1 : N}\right) \propto p\left(x_{1 : n} | \lambda_{1}\right) p\left(x_{n+1 : N} | \lambda_{2}\right) p\left(\lambda_{1}\right) p\left(\lambda_{2}\right) p(n)
$$

#### 条件付き分布の導出
これまでに説明したように，Gibbs Samplingでは確率変数の条件付き分布が必要となる．
条件付き分布の導出には，同時分布（Joint Distribution）と各変数間の依存関係から考えると良い．

パラメータ$n$，$\lambda_{1}$，$\lambda_{2}$の同時分布は次の通り：

$$
\begin{aligned}
p\left(x_{1 : n} | \lambda_{1}\right) p\left(x_{n+1 : N} | \lambda_{2}\right) p\left(\lambda_{1}\right) p\left(\lambda_{2}\right) p(n) \\
=\left(\prod_{i=1}^{n} p\left(x_{i} | \lambda_{1}\right)\right)\left(\prod_{i=n+1}^{N} p\left(x_{i} | \lambda_{2}\right)\right) p\left(\lambda_{1}\right) p\left(\lambda_{2}\right) p(n)
\end{aligned}
$$

上式の自然対数を取り，ポアソン・ガンマ分布の確率（密度）関数を代入した結果は次のようになる：

$$
\begin{equation}
\begin{array}{l}
\log p\left(x_{1 : n} | \lambda_{1}\right) +
\log p\left(x_{n+1 : N} | \lambda_{2}\right) +
\log p\left(\lambda_{1}\right) +
\log p\left(\lambda_{2}\right) +
\log p(n) \\
= \sum_{i=1}^{n}\left(x_{i} \log \lambda_{1}\right) - \lambda_{1} - \log\left(x_{i}!\right) \\ +
\sum_{i=n+1}^{N}\left(x_{i} \log \lambda_{2}\right) - \lambda_{2} - \log\left(x_{i}!\right) \\ +
\left(a-1\right) \log \lambda_{1} - b \lambda_{1} - \log \Gamma\left(a\right) + a \log b \\ +
\left(a-1\right) \log \lambda_{2} - b \lambda_{2} - \log \Gamma\left(a\right) + a \log b \\ -
\log N
\end{array}
\end{equation}
$$

これで各潜在変数の条件付き確率を獲得することができるようになった：**注目する確率変数を含む同時分布の項を集めれば良い**．

:::info
もともとの同時分布と自然対数をとった同時分布は上に示した通り．ある確率変数の条件付き分布は同時分布をそれ以外の確率変数の事前分布で除した形式で表されるが，対数を取ることでその操作は単なる引き算となる．その結果として，ある確率変数の条件付き分布はその確率変数を含む項を自然対数をとった同時分布から抽出するだけの作業と等価となるということ．

以下の変形を見れば，$\log p \left( b \right)$は，同時分布から$b$が独立して決める項を除去する操作と同じであることがわかる．

$$
\begin{aligned}
p\left(a | b\right) &= \frac{p\left(a, b\right)}{p\left(b\right)} \\
\log p\left(a | b\right) &= \log p\left(a, b\right) - \log p\left( b \right)
\end{aligned}
$$
:::

#### $\lambda_{1}$の条件付き事後分布
$\lambda_{1}$について注目すると，条件付き事後分布は次：
ただし$=^{+}$はある定数を加えれば等式が成り立つという意味．

$$
\begin{equation}
\begin{aligned} & \log p\left(\lambda_{1} | n, \lambda_{2}, x_{1 : N}\right) \\
=& \sum_{i=1}^{n}\left(x_{i} \log \lambda_{1}-\lambda_{1}\right)+(a-1) \log \lambda_{1}-b \lambda_{1} \\
=&\left(a+\sum_{i=1}^{n} x_{i}-1\right) \log \lambda_{1}-(n+b) \lambda_{1} \\
=& \log \operatorname{Gamma}\left(a+\sum_{i=1}^{n} x_{i}, n+b\right) \\
&+ \log\Gamma\left(a + \sum_{i=1}^{n}x_{i}\right) - \left(a + \sum_{i=1}^{n}x_{i}\right) \log \left(n + b\right) \\
=&^{+} \log \operatorname{Gamma}\left(a+\sum_{i=1}^{n} x_{i}, n+b\right) \end{aligned}
\end{equation}
$$

この結果は次のように読み解くことができる：
> 「$\lambda_{1}$についての条件付き事後分布は$a$→$a + \sum_{i=1}^{n}x_{i}$，$b$→$n + b$に置き換えたガンマ分布」

パラメータ$\lambda_{1}$について設定した事前分布と条件付き事後分布が同じになるような分布を共役事前分布という．

#### $\lambda_{2}$の条件付き事後分布
$\lambda_{2}$についても同様で，次のような条件付き事後分布が得られる：

$$
\begin{equation}
\begin{aligned} & \log p\left(\lambda_{2} | n, \lambda_{1}, x_{1 : N}\right) \\=& \sum_{i=n+1}^{N}\left(x_{i} \log \lambda_{2}-\lambda_{2}\right)+(a-1) \log \lambda_{2}-b \lambda_{2} \\=&^{+} \log \operatorname{Gamma}\left(a+\sum_{i=n+1}^{N} x_{i}, N-n+b\right) \end{aligned}
\end{equation}
$$

#### $n$の条件付き事後分布
$n$についての条件付き事後分布は次のようになる：

$$
\begin{equation}
\begin{array}{l}{\quad \log p\left(n | \lambda_{1}, \lambda_{2}, x_{1 : N}\right)} \\ {=\sum_{i=1}^{n}\left(x_{i} \log \lambda_{1}-\lambda_{1}-\log \left(x_{i} !\right)\right)+\sum_{i=n+1}^{N}\left(x_{i} \log \lambda_{2}-\lambda_{2}-\log \left(x_{i} !\right)\right)} \\ {=^{+}\left(\sum_{i=1}^{n} x_{i}\right) \log \lambda_{1}-n \lambda_{1}+\left(\sum_{i=n+1}^{N} x_{i}\right) \log \lambda_{2}-(N-n) \lambda_{2}}\end{array}
\end{equation}
$$

これは特に共役事前分布の形になっているわけではない．この場合には次の手順のような工夫を行う：
1. $n=1,\ldots,N$について上の式を計算
2. $\frac{\exp\left(n\right)}{\sum_{i}\exp\left(i\right)}$の演算で正規化
3. その結果を**多項分布**に当てはめてサンプリング → 新しい$n$とする

これらパラメータ更新とサンプリングをAlgorithmに基づいて繰り返すことでパラメータの事後分布を得られる．

### Sample Code

from scipy.stats import uniform, gamma, poisson
from matplotlib import pyplot as plt
import numpy as np

np.random.seed(123456789)

N = 50
a = 2  # param of Gamma Dist
b = 1  # param of Gamma Dist

# 変化点
n = int(round(uniform.rvs() * N))
print(f'変化点: {n}')

# 各期間での頻度の平均値
# 1./bとしているのはrandomの引数の仕様がそうなっているため
lambda1 = gamma.rvs(a, scale=1./b)
lambda2 = gamma.rvs(a, scale=1./b)

lambda_series = np.ones(N)
lambda_series[:n] = lambda1
lambda_series[n:] = lambda2

# カウントデータ
x = poisson.rvs(lambda_series)

# Store the samples
E = 5200
BURN_IN = 200

chain_n = np.zeros(E - BURN_IN)
chain_lambda1 = np.zeros(E - BURN_IN)
chain_lambda2 = np.zeros(E - BURN_IN)

# Gibbs Sampling
for e in range(E):
    if e % 1000 == 999:
        print(f"At iteration {e + 1}")
    # sample lambdas
    lambda1 = gamma.rvs(a + np.sum(x[:n]), scale=1./(n + b))
    lambda2 = gamma.rvs(a + np.sum(x[n:]), scale=1./(N - n + b))
    
    # sample n
    mult_n = np.zeros(N)
    for i in range(N):
        mult_n[i] = np.sum(x[:i]) * np.log(lambda1) - i * lambda1 + np.sum(x[i:]) * np.log(lambda2) - (N - i) * lambda2
        
    # 多項分布のパラメータを作ってサンプリング
    mult_n = np.exp(mult_n - np.max(mult_n))  # 多項分布のパラメータにする
    n = np.where(np.random.multinomial(1, mult_n / mult_n.sum(), size=1)==1)[1][0] # サンプリング結果
    
    if e >= BURN_IN:
        idx = e - BURN_IN
        chain_n[idx] = lambda1
        chain_lambda1[idx] = lambda2

# Plotting
# plt.figure(figsize=(32, 18))
f, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, 1)

ax1.stem(range(N), x, linefmt='b-', markerfmt='bo')
ax1.plot(range(N), lambda_series, 'r--')
ax1.set_ylabel('Counts')

ax2.plot(chain_lambda1,'b',chain_lambda2,'g')
ax2.set_ylabel('$\lambda$')

ax3.hist(chain_lambda1,20)
ax3.set_xlabel('$\lambda_1$')
ax3.set_xlim([0,12])

ax4.hist(chain_lambda2,20,color='g')
ax4.set_xlim([0,12])
ax4.set_xlabel('$\lambda_2$')

ax5.hist(chain_n,50)
ax5.set_xlabel('n')
ax5.set_xlim([1,50])

plt.show()

## Reference
- [Bayesian Inference: Gibbs Sampling](http://www.mit.edu/~ilkery/papers/GibbsSampling.pdf)
- [Qiita: 変化点検出入門](https://qiita.com/nozma/items/a27cf0d327014ad1ae0f)

## Appendix
### ベイズ推論のワークフロー
1. 応答変数のモデリング（確率分布の当てはめ）<br />→ $p\left(\mathcal{D} | \theta\right)$
2. 確率変数間の依存関係を整理し，生成過程をグラフィカルモデルで記述する<br />→ $p\left(\theta\right) = p\left(\theta | \alpha\right) p\left(\alpha\right)$
3. 確率変数の同時分布を依存関係に沿って分解し，（考えやすいように）自然対数を取る<br />→ $\log p\left(\theta, \mathcal{D}\right)$
4. 事後分布を同時分布と事前分布から求める<br />→ $\log p\left(\theta | \mathcal{D}\right) = \left[\log p\left(\theta | \mathcal{D}\right) - \log p\left(\mathcal{D}\right)\right] \propto \log p\left(\theta, \mathcal{D}\right)$
5. 4の結果から新たにサンプリングし更新する確率変数を含む項だけを抽出，その結果は条件付き分布となる<br />→ $\log p\left(\theta, \mathcal{D}\right) =^{+} \theta_{i} \text{related terms}$
6. 条件付き分布を利用してサンプリングと更新を行う