「基礎からのベイズ統計学」を読んで突っかかった点

# 目的
「基礎からのベイズ統計学」を読んでます

https://www.asakura.co.jp/detail.php?book_code=12212

4章MH法で突っかかってしまったので、悩んだ点と解決をまとめていきます。

# 離散マルコフ連鎖の定常条件
## マルコフ連鎖について
$X_1,\cdots,X_n$と確率変数が順々更新されていく状況を考えます。
確率変数がとり得る状態は$a,b$のみであるとします。

ここで、更新則は

```math
\left[\begin{array}{c}
P\left(X_2=a\right) \\
P\left(X_2=b\right)
\end{array}\right]
=
\left[\begin{array}{c}
P\left(X_2=a\right|X_1 = a) &P\left(X_2=a\right|X_1 = b)\\
P\left(X_2=b\right|X_1 = a) &P\left(X_2=b\right|X_1 = b)
\end{array}\right]

\left[\begin{array}{c}
P\left(X_1=a\right) \\
P\left(X_1=b\right)
\end{array}\right]
```
になります。

ここで

```math
\left[\begin{array}{c}
P\left(X_2=a\right|X_1 = a) &P\left(X_2=a\right|X_1 = b)\\
P\left(X_2=b\right|X_1 = a) &P\left(X_2=b\right|X_1 = b)
\end{array}\right]
```
の列和は1になります。

```math
\left[\begin{array}{c}
P\left(X_2=a\right|X_1 = a) &P\left(X_2=a\right|X_1 = b)\\
+&+\\
P\left(X_2=b\right|X_1 = a) &P\left(X_2=b\right|X_1 = b)
\end{array}\right] \\
=\\
\left[\begin{array}{c}
1 & 1
\end{array}\right] \\
```

ここで遷移行列を$P$と書けば

```math
\left[\begin{array}{c}
1 & 1
\end{array}\right] \\
P
=
\left[\begin{array}{c}
1 & 1
\end{array}\right] \\
```
ということなので、$P$は固有値1を持ちます。$P^T$の固有ベクトルは

```math
\left[\begin{array}{c}
1 \\
\vdots\\
1\\
\end{array}\right]
```
を持つこともわかります。

```math
\left[\begin{array}{c}
P\left(X_2=a\right|X_1 = a) &P\left(X_2=a\right|X_1 = b)\\
P\left(X_2=b\right|X_1 = a) &P\left(X_2=b\right|X_1 = b)
\end{array}\right]
```
の行和は

```math
\left[\begin{array}{c}
P\left(X_2=a\right|X_1 = a) + P\left(X_2=a\right|X_1 = b)\\
P\left(X_2=b\right|X_1 = a) + P\left(X_2=b\right|X_1 = b)
\end{array}\right]

=

\left[\begin{array}{c}
P\left(X_2=a\right)\\
P\left(X_2=b\right)
\end{array}\right]
```
となるので、$P$は固有ベクトルとして

```math
\left[\begin{array}{c}
1 \\
\vdots\\
1\\
\end{array}\right]
```
を持つとは限らないこともわかります。

## マルコフ連鎖の定常条件
遷移則は

```math
PX_{n-1} = X_n
```
であります。

```math
X_n =\lambda_1 v_1 + \cdots + \lambda_n v_n
```
と固有ベクトルの合成で作られると考えてみましょう。

```math
P^m X_n = P^m\lambda_1 v_1 + \cdots + P^m \lambda_n v_n
```
となって、最大固有値支配と考えられます。つまり

```math
P^m X_n = \lambda_1^m v_1 + \cdots + \lambda_n^m v_n
```
となって、固有値が1未満なら0に収束し、1を超えるなら発散していきます。1なら定常状態になります。
マルコフ連鎖の遷移行列は固有値１を必ず持つので0に収束することはありえません。
発散することも、確率ベクトルであることからありえません。
よって、**定常状態に収束したときの確率ベクトルは固有値1の固有ベクトルです。**
固有値1に複数の固有ベクトルがある可能性もあります（追記：ペロン・フロベニウスの定理からありえないみたいです）。

## 固有値が1を超えない証明
確率の話で条件付き確率で更新していくのがマルコフ連鎖です。確率ベクトルが発散するという状態はありえないと言えます。ただ、定性的ではなく行列の性質から導いてみましょう。

https://math-note.xyz/linear-algebra/perron-frobenius-theorem/

ペロン＝フロベニウスの定理で示されるみたいです。
読むのはだいぶきついのでやめました。

## マルコフ連鎖が定常条件

```math
\left[\begin{array}{c}
P\left(X_{n+1}=a\right) \\
P\left(X_{n+1}=b\right)
\end{array}\right]
=
\left[\begin{array}{c}
P\left(X_{n+1}=a\right|X_n = a) &P\left(X_{n+1}=a\right|X_n = b)\\
P\left(X_{n+1}=b\right|X_n = a) &P\left(X_{n+1}=b\right|X_n = b)
\end{array}\right]

\left[\begin{array}{c}
P\left(X_n=a\right) \\
P\left(X_n=b\right)
\end{array}\right]
```
の更新則で変化しないということなので、


```math
\forall a\\
について
P(X_n = a) = \sum_{b \in N} P(X_{n+1} = a | X_n=b) P(X_n =b)
が成立する。
```
これがマルコフ連鎖の定常条件です。

https://manabitimes.jp/math/2757


# 連続のマルコフ連鎖とは
次に連続確率密度関数を持つ1確率変数のマルコフ連鎖を考えます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/fae1ce56-2269-3553-a4d5-746f2eaa8caf.png)
確率変数$X_1$は確率密度関数$f(X_1)$に従うとします。
$X_1$が特定の値を取った後、確率変数$X_2$は確率密度関数$f(X_2|X_1=a)$に従うと考えます。
これが、連続確率密度関数を持つ１つの確率変数のマルコフ連鎖です。

もう少し具体的に考えてみましょう。

## 正規分布のマルコフ連鎖

```math
X_1 \sim f(X_1) = N(0,1)\\
```

```math
X_2 \sim f(X_2|X_1) = N(x_1,1)
```
に従うと考えます。
$X_2$の分布は一つ前の$X_1$の結果によって変わります。
一つ前の値を中心に正規分布するという感じです。
ランダムウォークみたいですね。

漸化式で書けば

```math
X_n \sim f(X_n|X_{n-1}) =N(x_{n-1},1)
```
ここで、周辺確率の求め方で行けば

```math
f(X_2) = \int f(X_2 | X_1)f(X_1)dX_1
```
になります。$f(X_1)$が事前分布で、$f(X_2 | X_1)$が遷移核といえます。

離散確率を持つ確率変数の遷移行列Pのようなものは定義できない（多変数にしたら定義できるかも）ですが、条件付き確率（密度）が遷移核・遷移行列であることは同じです。

## 実際に遷移則を求めてみる
先程の正規分布の例で遷移則を求めてみましょう。

```math
f(X_2) = \int f(X_2 | X_1)f(X_1)dX_1 = \int N(x_1,1) N(0,1) dx_1 = \int \frac{1}{\sqrt{2 \pi}} exp(- \frac{(x_2-x_1)^2 + x_2^2}{2})dx_1
```
$f(X_2)$は$X_1$が未知の場合の$X_2$の分布になります。未知というと、観測者が知っているかどうかみたいな話になってしまいますが、観測した時点で特定の分布の存在が否定されるだけと考えてください。
右辺の積分を計算すれば新たな事前分布が出来上がります。今度は$f(X_3)$を求めてみましょう・・・というのが連続確率変数のマルコフ連鎖になります。積分計算が入るので実際に行うのは難しいです。
積分の過程で一つ前の変数が消えることも重要です。
この計算過程が遷移則です。

## 定常条件

```math
\forall a, \forall b に対して
f(a | b)f(b) = f(b|a)f(a)
```
が成立していると仮定すると、マルコフ連鎖の遷移によって発生する事前分布と事後分布が一致します。

```math
f(X_n) =f(X_{n-1})
```

### 証明

確率密度関数を$f$として記述してきましたが、事後・事前確率、条件付き確率で同じ文字が並ぶため分かりにくくなります。そこで、次のように名前を変えておきます。

```math
f(X_n) = f_n(X_n)\\
```

```math
f(X_{n-1}) = f_{n-1}(X_{n-1})\\
```
```math
f(X_n|X_{n-1}) = g(X_n|X_{n-1})\\
```
名前を変えただけです。
証明をしていきます。

```math
f(X_n) = \int f(X_n | X_{n-1})f(X_{n-1})dX_{n-1}
```
と周辺確率の計算からわかっています。この式は$X_n = a,X_{n-1} = b$と特定の値を値を確率変数が取った場合も成立する式です。
先程の名前の変更も併せて書き換えると

```math
f_n(a) = \int g(a | b)f_{n-1}(b)db
```
となります。定常条件の式を再度上げると次になります。

```math
\forall a, \forall b に対して
f(a | b)f(b) = f(b|a)f(a)
```

先程の名前の書き換えを適用すると

```math
\forall a, \forall b に対して
g(a | b)f_{n-1}(b) = g(b|a)f_{n-1}(a)
```

これを、使って書き換えます！

```math
f_n(a) =  \int g(b|a)f_{n-1}(a)db = f_{n-1}(a) \int g(b|a)db = f_{n-1}(a)
```
以上より、

```math
f(X_n) = f(X_{n-1})
```
を示すことができました。

なぜ名前の書き換えをしたかというと、書き換えずに定常条件を書くと

```math
f(X_n |X_{n-1})f(X_{n-1}) = f(X_{n-1}|X_n)f(X_n)
```
になりますが、右辺の$f(X_n)$と左辺の$f(X_{n-1})$はともに$X_{n-1}$の確率密度関数です。$f$の文字が同じ為勘違いしてしまいます（原因追求に3日かかった）。

# MH法を実践

MH法の目的、**指定した分布従う確率変数から実際に値を複数取得すること。**

## 目的の目的
任意の分布から平均値や分散を計算するのは積分計算が入るため、実行不可能な場合がある。


## 問題設定

マルコフ連鎖により定常分布になっているので

```math
(事前分布) = (事後分布)
```

問題設定として、

```math
(事後分布)のカーネルは既知
```
であるとしています。カーネルのみでいいのは、事後分布と事前分布の割り算をするため。

```math
f(X_2) = \int f(X_2 | X_1)f(X_1)dX_1
```
で遷移していき、定常になったとき

```math
f(X) = \int f(X_2 | X)f(X)dX
```
$f(X)$は事前にわかっている状況で、定常条件より

```math
\forall a, \forall b に対して
f(a | b)f(b) = f(b|a)f(a)
```
が成立している必要がある。

## MH法
遷移確率$f(a|b)$はわからないので、とりあえず計算者側が任意に設定する。
任意に設定する以上、定常条件を満たすとは限らない。

```math
f(a | b)f(b) > f(b|a)f(a)
```
を判定する。(>,<,=全てありえる)

```math
r= \frac{f(b|a)f(a)} {f(a | b)f(b) }
```
として、

```math
r f(a | b)f(b) =  f(b|a)f(a)
```
$r$をかければ、等号成立します。定常条件が成立しているということです。


## 実行

1. 任意に$b$を設定（乱数的に生成）、
1. 遷移確率分布$f(|b)$に従い、$a$を生成。$f(a|b)$を計算。
1. $f(b|a)$を計算
1. $f(a)$を計算
1. $f(b)$を計算
1. $r$を計算

定常条件

```math
r f(a | b)f(b) =  f(b|a)f(a)
```
を

```math
(r f(a | b)) *f(b) =  f(b|a)f(a)
```
として

```math
f^\prime (a|b) = r f(a | b)
```
と考えてみましょう。元々、推移確率分布はわからないので任意で指定していました。ので、これを新たな推移確率と見てみましょう。

```math
bを元に生成したaの採用確率を r f(a | b)にする。
```
ということです。仮に$r>1$なら100%$a$を採用します。
採用する場合は、$a$を新たな$b$として、前述の計算を繰り返します。
以上を無限に繰り返せば、


```math
\forall a, \forall b に対して
f(a | b)f(b) = f(b|a)f(a)
```
が成立していると考えられます。

## 補足：rの計算

```math
r= \frac{f(b|a)f(a)} {f(a | b)f(b) }
```
であったので、

```math
r= \frac{f(b|a)f(a)_{カーネル}} {f(a | b)f(b)_{カーネル} }
```
で計算できます。つまり、事後分布はカーネルがあれば計算できます。


## Python
実際にやってみます
先に言えば、

https://qiita.com/Fumio-eisan/items/6e8bad9977e945ddb2fc

を参考にしたほうがいいです。

推移確率分布として正規分布を使います。
選定条件としては、分布に従う乱数の発生ができること・確率を計算できることが必要です。

```math
f(a|b) = N(平均b,分散1)
```
とします。

正規分布に従う乱数を取得

https://biotech-lab.org/articles/13652

```python
import random
import numpy as np
import matplotlib.pyplot as plt
import math

b=0

import random
a = random.gauss(b,1)
```

推移確率を計算（カーネル値）

```python
f_ab = math.exp(-1*(a-b)**2/2)
f_ba = math.exp(-1*(b-a)**2/2)
```

求めたい事後分布をガンマ分布として、カーネルを

```math
f(x) = e^{-13x} x^{10}
```
とします。

生成確率を計算(値がでかすぎるとエラーになったため分岐)

```python
def f(x):
    if abs(x) > 10:
        f_x = 1e+10
    else:
        f_x = math.exp(-13*x)*(x**10)        
    return f_x
```

$r$を計算。0割が発生する可能性があるため、条件分岐させておきます。

```python
if f_b > 0.00000000001 and f_ab >0.00000000001:
    r = f_ba*f_a/f_ab/f_b
else:
    r = 1
```

$a$を採用するかどうか判定＋ループにする。

```python
rnd = np.empty([0]) #1次元から配列生成:保存用
b = 0
for i in range(100000):
    
    a = random.gauss(b,1)
    
    f_ab = math.exp(-1*(a-b)**2/2)
    f_ba = math.exp(-1*(b-a)**2/2)
    
    f_a = f(a)
    f_b = f(b)
    
    if f_b > 0.00000000001 and f_ab >0.00000000001 :
        r = f_ba*f_a/f_ab/f_b
    else:
        r = 1
    
    print(r)
    
    if r > 1 :
        b = a
        rnd = np.append(rnd,b)
    else:
        hoge = random.uniform(0,1)
        if hoge < r:
            b = a
            rnd = np.append(rnd,b)
```


ヒストグラムを表示(密度分布に変換)

```python
# ヒストグラム
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(rnd, bins=100, histtype='barstacked', ec='black')
plt.show()
```

最終的に

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/90ddac43-442b-664b-1d10-4e87a5252c3b.png)

```python
import random
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.stats as st


def f(x):
    """
    if abs(x) > 20:
        f_x = 1e+10
    else:
        f_x = math.exp(-13*x)*(x**10)
        
    if f_x < 0:
        f_x = 0
    """
    
    f_x = math.exp(-13*x)*(x**10)
    
    return f_x



rnd = np.empty([0]) #1次元から配列生成:保存用
b = 2
for i in range(1*10**6):
    
    while True:
        a = random.gauss(b,0.5)
        if a > 0:
            break
    
    
    f_ab = math.exp(-1*(a-b)**2/0.5)
    f_ba = math.exp(-1*(b-a)**2/0.5)
    
    f_a = f(a)
    f_b = f(b)
    
    r = f_ba*f_a/f_ab/f_b
    
    """
    if abs(f_b*f_ab) > 1e-6 :
        r = f_ba*f_a/f_ab/f_b
    else:
        r = 1
    """
    
    #print(r)
    
    if r > 1 :
        b = a
        rnd = np.append(rnd,b)
    else:
        hoge = random.uniform(0,1)
        if hoge < r:
            b = a
            rnd = np.append(rnd,b)
            
                
            
        
    
# ヒストグラム
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(rnd, bins=20, histtype='barstacked', ec='black',density=True)

xx = np.linspace(0, 2.5,501)
ax.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))

plt.show()    

print("f_ab "+str(f_ab))
print("f_b "+str(f_b))
print("r  "+str(r))
print("---size---")
print(rnd.size)
```
このグラフが得られるかどうかは、かなり初期値に依存します。


# ランダムウォークMH法

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/2eb8b0a0-5f1c-8ec9-7df2-911cd656a5d3.png)

ループを増やした結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/116f1e30-cf82-f780-2af7-b103b1ef0d38.png)

先端があいませんね

```python

import random
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.stats as st


def f(x):
    """
    if abs(x) > 20:
        f_x = 1e+10
    else:
        f_x = math.exp(-13*x)*(x**10)
        
    if f_x < 0:
        f_x = 0
    """
    
    f_x = math.exp(-13*x)*(x**10)
    
    return f_x



rnd = np.empty([0]) #1次元から配列生成:保存用
b = 2
for i in range(1*10**5):
    
    while True:
        a = b + random.gauss(0,0.5)
        if a>0:
            break
    
    
    #f_ab = math.exp(-1*(a-b)**2/0.5)
    #f_ba = f_ab
    
    f_a = f(a)
    f_b = f(b)
    
    #r = f_ba*f_a/f_ab/f_b
    r = f_a/f_b
    
    """
    if abs(f_b*f_ab) > 1e-6 :
        r = f_ba*f_a/f_ab/f_b
    else:
        r = 1
    """
    
    #print(r)
    
    if r > 1 :
        b = a
        rnd = np.append(rnd,b)
    else:
        hoge = random.uniform(0,1)
        if hoge < r:
            b = a
            rnd = np.append(rnd,b)
            
                
            
        
    
# ヒストグラム
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(rnd, bins=100, histtype='barstacked', ec='black',density=True)

xx = np.linspace(0, 2.5,501)
ax.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))

plt.show()    

#print("f_ab "+str(f_ab))
print("f_b "+str(f_b))
print("r  "+str(r))
print("---size---")
print(rnd.size)
```

# 最後に
まだ読みきれていないので、発見があり次第まとめたいと思います。