動機

MCMCで計算をすると，カルマンフィルタを実装したときは実感できたベイズ更新の部分がよくわからないという僕の気持ちの問題がありました。
そこで，須山さんの「ベイズ推論による機械学習入門」

機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)

作者: 須山敦志,杉山将
出版社/メーカー: 講談社
発売日: 2017/10/21
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

の3章で，学習と推論を解析的に導出していたため，自分で計算・実装を行い理解を深めたいというのがこの記事を書く動機です。

線形回帰の式

１次関数での線形回帰

想像しやすいように，１次関数での線形回帰を行います。式は $y = a x + b+ \varepsilon_r$ です。
$\varepsilon_r$ はノイズで，回帰分析の事後分布として求めたい値は， $a, b$ となります。
ノイズは簡単のために正規分布に従うと仮定して，
$\varepsilon_r \sim N(\varepsilon_r | 0, \lambda^{-1})$
と定義します。テキストに習って $\lambda$ は既知であるとします。
ある一つの観測値 $y_k$ とある一つの入力値 $x_k$ が得られたときに，事後分布 $p(a| y_k, x_k), p(b|y_k, x_k)$ が導出できれば，学習と予測ができることになります。

まず，簡単のために，bは既知であるとして，aの事後分布を解析的に求めてみます。
ベイズの定理より，求めたい事後分布 $p(a| y_k, x_k)$ は，
$p(a| y_k, x_k) = \frac{p(y_k|x_k, a) p(a)}{p(y_k, x_k)}$
と書けます。
このときに，事前分布 $p(a)$ は正規分布していると仮定し
$p(a) = N(a | m_a, \Lambda_a^{-1})= \frac{1}{\sqrt{2\pi\Lambda_a^{-1}}}\exp\big( -\frac{1}{2}\Lambda_a(a-m_a)^{2} \big)$
と書けるとします。ここで， $m_a, \Lambda_a$ は事前分布の平均と分散で固定値で与えるハイパーパラメータです。

さて，事前分布 $p(a)$ と入力値 $x_k$ が与えられた元で，観測値 $y_k$ が得られる度合いを表す尤度を考えます。
先程，ノイズが正規分布に従うとしたため，
$p(y_k|x_k, a) = N(y_k | a * x_k +b, \lambda^{-1} )$
と書くことができます。ここで，この正規分布を真面目に書いてみると
$N(y_k | a * x_k+b, \lambda^{-1}) = \frac{1}{\sqrt{2\pi\lambda^{-1}}}\exp\big( -\frac{1}{2}\lambda(y_k-a*x_k-b)^{2} \big)$
となります。
今は $a$ の事後分布だけが気になっているので，計算を簡単にするために事前分布，尤度をそれぞれの対数( $\log$ )を取って， $a$ について整理していきます。

事前分布

$\log p(a) = \log \frac{1}{\sqrt{2\pi\Lambda_a^{-1}}}\exp\big( -\frac{1}{2}\Lambda_a(a-m_a)^{2} \big)$
$= -\frac{1}{2} \log (2\pi\Lambda_a^{-1} ) -\frac{1}{2}\Lambda_a(a-m_a)^{2}$

第一項は $a$ に関係のない項なので，第二項を展開して， $a$ に関する項を昇順に並べていきます。
$-\frac{1}{2}\Lambda_a(a-m_a)^{2} = -\frac{1}{2} \Lambda_a a^{2} + \Lambda_a m_a a -\frac{1}{2} \Lambda_a m_a^{2}$ (1)
となります。

尤度

次に，先程定義した尤度も同様に計算することで，
$-\frac{1}{2} \lambda_a x_k ^{2} a^{2} + \lambda (y_k -b) x_k a -\frac{1}{2} \lambda (y_k -b )^{2}$ (2)
と書けます。

事前分布と尤度の積

ベイズの定理の分子である $p(y_k|x_k, a) p(a)$ は事前分布と尤度の積ですので，対数を取ったあとは和になります。
なので，(1)と(2)の和を計算して， $a$ について昇順に並べることで，事後分布の $a$ についてはどのような形になるのか想像することができます。和を取ると $-\frac{1}{2} ( \Lambda_a + \lambda_a x_k ^{2} ) a^{2} +( \Lambda_a m_a + \lambda (y_k -b) x_k )a -\frac{1}{2} ( \Lambda_a m_a^{2} + \lambda (y_k -b )^{2})$ (3)
と書けます。
この形は事前分布や尤度で仮定した正規分布と全く同じ形になっていることから，事後分布も正規分布であると考えることができます。*1
分母は，周辺尤度を求めれば，正規化されていると思います。
このあたりの数式は天下り的に計算はしませんが，任意の事後分布の正規化の数式などは例えば，

ノンパラメトリックベイズ点過程と統計的機械学習の数理 (機械学習プロフェッショナルシリーズ)

作者: 佐藤一誠
出版社/メーカー: 講談社
発売日: 2016/04/20
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

こちらの本などに書かれているように思えます。
さて，事後分布も正規分布であるとして，一点を観測したあとの傾き $a$ の事後分布は
$p(a| y_k, x_k) = N(a| ( \Lambda_a + \lambda_a x_k ^{2} ), \frac{\Lambda_a m_a + \lambda (y_k -b) x_k}{( \Lambda_a + \lambda_a x_k ^{2} )})$ (4)
として計算できます。*2

今，観測点が一点だけの場合を考えていましたが，複数観測する場合には，式(3)は更に一般化することができて， $-\frac{1}{2} ( \Lambda_a +\sum_k \lambda_a x_k ^{2} ) a^{2} +( \Lambda_a m_a + \sum_k \lambda (y_k -b) x_k )a -\frac{1}{2} ( \Lambda_a m_a^{2} + \sum_k \lambda (y_k -b )^{2})$ (3)
と，書けます。このとき，(4)は
$p(a| y_k, x_k) = N(a| \frac{\Lambda_a m_a + \sum_k \lambda (y_k -b) x_k}{( \Lambda_a + \sum_k \lambda_a x_k ^{2} )}, ( \Lambda_a + \sum_k \lambda_a x_k ^{2} ))$ (5)
となりますので，この式を使うことで，一点一点を更新していく様子がみることができます。