本当に正規分布同士の積で分布は更新されるのか

ふと，確率ロボティクスを見ていたとき，今更ながら事前分布と尤度の積は本当に事後分布の形になるのかが気になった．

今まで当たり前のように「更新される」と考えていたが，はたして本当にそうなのだろうか？
というか，正規分布の積は本当に積分布になるのか？
数学的証明はできるけど，数学にそこまで強くないので，直観で理解したい．

なので，今回はその当たり前のことをプログラムにして確認した．

具体的には
正規分布のグラフを2つ作って，その2つを掛け算して更新したグラフと
平均と分散を最初に計算して求めたグラフが一致するのかを調べた．

グラフの作成

正規分布の式は

\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}

で表されるため

import math
import numpy as np
import matplotlib.pyplot as plt

def gaussian(mu, sigma, x):
  return np.exp(-(x-mu)**2/(2*sigma**2)) / (math.sqrt(2*math.pi)*sigma)

# 正規分布1
mu1 = 4
sigma1 = 3

# 正規分布2
mu2 = -3
sigma2 = 2

# グラフのx軸の幅
step = 100000
x_lim = 10.0

n = np.linspace(-x_lim, x_lim, step)
# グラフの1メモリの横幅
interval = x_lim*2 / step


# 正規分布をつくる
prior_distribution1 = gaussian(mu1, sigma1, n) # 正規分布1
prior_distribution2= gaussian(mu2, sigma2, n) # 正規分布2
posterior_distribution = gaussian(mu1, sigma1, n) * gaussian(mu2, sigma2, n) # 正規分布1と正規分布2の積


# 正規化定数（面積を1にするための定数）を求める．1.0÷面積で求まる
normalization = 1.0 / (np.array(posterior_distribution).sum() * interval)
# 正規化
posterior_distribution = np.array(posterior_distribution) * normalization

# グラフに表示
plt.plot(n, prior_distribution1, label="prior_distribution1")
plt.plot(n, prior_distribution2, label="prior_distribution2")
plt.plot(n, posterior_distribution, label="posterior_distribution")
plt.xlim([-x_lim, x_lim])
plt.legend()
plt.show()

結果として

二つの正規分布の積が得られていることを確認した．

グラフが正しいか確認

計算で平均と分散を求める

始めに2つの正規分布の平均と分散から，事後分布の平均と分散を導いてみる．

「道具としてのベイズ統計」p132より
正規分布1の平均と分散を$\mu1$と$\sigma1$
正規分布2の平均と分散を$\mu2$と$\sigma2$とすると，

平均は

\mu = \frac{\frac{\mu1}{{\sigma1}^{2}} + \frac{\mu2}{{\sigma2}^{2}}} {\frac{1}{{\sigma1}^{2}} + \frac{1}{{\sigma2}^{2}}}

分散は

{\sigma}^{2} = \frac{1} {\frac{1}{{\sigma1}^{2}} + \frac{1}{{\sigma2}^{2}}}

より求められる．

この式をpythonに起こして計算すると

# 平均
mu3 = ((mu2/sigma2**2) + (mu1/sigma1**2)) / ((1/sigma2**2) + (1/sigma1**2))
# 分散
sigma3 = 1.0 / ( (1.0/sigma1**2.0) + (1.0/sigma2**2.0) )
sigma3 = sigma3**0.5

print("計算で求めた平均値")
print(mu3)
print("計算で求めた分散")
print(sigma3)

計算で求めた平均値
-0.8461538461538463
計算で求めた分散
1.6641005886756874

グラフから平均と分散を求める

グラフの最大値より平均を求めると

print("グラフから求めた平均")
print(n[np.argmax(posterior_distribution)])

グラフから求めた平均
-0.8461084610846097

次に，グラフから分散を求める

分散はグラフ全体の面積の68.3%のとき，x軸が$2\sigma$となるから，$\mu$から面積の34.15%までのx軸の値を調べることで，$\sigma$を求める（精度に問題あり）．

area_sum = 0
for i in range(int(step/2)):
  area_sum += posterior_distribution[mu3_num + i] * interval
  if area_sum >= 0.3415: # 68.3%
    print("グラフから求めた分散")
    print(i* interval)
    break

グラフから求めた分散
1.665

グラフと計算での平均と分散を比較

グラフから求めた平均
-0.8461084610846097
計算で求めた平均値
-0.8461538461538463

グラフから求めた分散
1.665
計算で求めた分散
1.6641005886756874

あまり誤差が無いことが分かる（ここら辺は変数stepやx_limの値を大きくすればもっと正確になる）．

グラフの重ね合わせ

さきほどプログラムでつくった積のグラフと，計算で求めたグラフを合わせると

# 正規分布をつくる
posterior_distribution2 = gaussian(mu3, sigma3, n)

# グラフに表示
plt.plot(n, posterior_distribution, label="posterior_distribution", color='r')
plt.plot(n, posterior_distribution2, '--', label="posterior_distribution2", color='b')
plt.xlim([-x_lim, x_lim])
plt.legend()
plt.show()

ぴったり！

ということで，結果として
2つを掛け算して更新したグラフと，計算により求めたグラフが一致して
正規分布同士の積で正しく分布が更新されることが分かった．