---
disqus: ierosodin
---

# From LSE to Fully bayesian

> Organization contact [name= [ierosodin]([e-mail])]
###### tags: `machine learning` `學習筆記`

==[Back to Catalog](https://hackmd.io/@ierosodin/Machine_Learning)==

## LSE

回顧前面的章節，Linear regression 就是希望找出一條線來 fit 我們的 data，所以我們去找 least square error 的地方，那現在我們有了 bayes，也認識了高斯分佈，是否能將機率的概念用在這裡呢？也就是從前面的 ==frequentist 轉換成 bayesian==。

## MLE

沒錯，其實就是用前面所提到的 MLE，我們希望找到一條線，對於每一個 x，我們都能在線上產生一個高斯分佈，並利用計算 maximum likelihood，找出最有可能的那調線，如下圖。


其中，${\vec{w}}$ 跟前面一樣，是我們希望得到的那條線的參數：
${y\ =\ w_0\ +\ w_1x\ +\ w_2x^2\ +\ ...}$

然後我們一樣有一個 D 維的 design matrix，D 代表我們的資料量：
${y\ =\ Xw\ =\ \left[\begin{array}{ccccc}
1 & x_1 & x_1^2 & ... & x_1^k \\
1 & x_2 & x_2^2 & ... & x_2^k \\
1 & x_3 & x_3^2 & ... & x_3^k \\
... & ... & ... & ... & ... \\
1 & x_D & x_D^2 & ... & x_D^k
\end{array}\right]\left[\begin{array}{c}
w_0 \\
w_1 \\
w_2 \\
... \\
w_k
\end{array}\right]}$

每一個 x 上都會有一個高斯分佈，有其各自的 mean 與 variance，那如何判斷一個好的 ${w}$ 呢？就是找其 MLE：
${P(D|w)\ =\ P(Y|X,w)\ =\ \prod_i\frac{1}{\sqrt{2\pi}\sigma_i}e^{\frac{-1}{2\sigma_i^2}(y_i-x_iw)^2}\ \propto\ e^{\sum_i\frac{-1}{2\sigma_i^2}(y_i-x_iw)^2}}$

和前面一樣，連乘找極值就是取 ${log}$：
${maxP(D|w)\ =\ \sum_i\frac{-1}{2\sigma_i^2}(y_i-x_iw)^2\ \propto\ \sum_i(y_i-x_iw)^2}$

這裡可以發現，跟前面的 ${||A\vec{x}\ -\ \vec{b}||^2}$ 根本一樣。

## MAP

我們都知道，如果只使用 LSE，可能會產生 overfitting 的問題，因此當時利用增加一項 ${\lambda ||w||}$，來對係數做懲罰。

到了 Bayesian 一樣有類似的問題，今天如果只用 MLE，那我對 ${w}$ 的 knowledge 是 0，也就是今天有麼資料我就算甚模，找出這筆資料最適合的 ${w}$，至於==其他 ${w}$ 我都認為是不可能==，同樣也會造成 overfitting。

那在 MLE 要怎麼對係數做懲罰呢？其實就是用到 Bayes，也就是今天我有過去累積下來的 prior，告訴我哪些 ${w}$ 比較有可能是我要的，至於其他的 ${w}$，我==仍給予一定的出現機率==。換個方式說，我給一個 ${P(w)}$ 為高斯分佈，告訴我在過去的 knowledge 中，甚模 ${w}$ 是最容易出現的，再搭配新進來的 data，對 prior 做更新，得到新的 posterior：
${P(w|D)\ =\ \frac{P(D|w)P(w)}{marginal}}$

有了 prior 與 likelihood，我們一樣對它做微分取極值：
${max\ =\ log(P(D|w)P(w))\ =\ log(e^{\sum_i\frac{-1}{2\sigma_i^2}(y_i-x_iw)^2}e^{\frac{-1}{2}w^TbIw})}$

其中，${b}$ 為 prior 的 variance 的倒數，並且假設 mean 為零 (初始)。這裡再假設 ${a}$ 為 likelihood 中，variance 的倒數。

${max\ =\ log(e^{\frac{-a}{2}(y_i-x_iw)^2+\frac{-1}{2}w^TbIw})\ =\ \sum_i\frac{-a}{2}(y_i-x_iw)^2+\frac{-b}{2}w^Tw\ =\ \frac{-a}{2}(||Xw-y||^2\ +\ \frac{b}{a}w^Tw)}$

後面的 ${\frac{b}{a}}$ 其實就是所謂的懲罰項，與 rLSE 不同的地方在於，==在 rLSE 中，懲罰項都是相對於原點==，也就是希望係數 ${w}$ 不要太大，但==實際的 ${w}$ 位置可能與 0 相聚很遠==，因此 ${\lambda}$ 怎麼取會影響到結果的正確性，但這裡我是利用 likelihood 來更新 prior，因此==即使一開始假設的 mean 與 variance 是錯的，也能透過多次迭代後，往正確的方向移動==，而不會因為一開始假設錯誤就永遠錯誤。

> MAP 其實就是用到前面的 Bayesian classifier，只是前面提到的例子中，likelihood 是離散的，而這裡對於每一個 x，都會有一個連續的 gaussian distribution，同樣可以利用 prior 與 likelihood 去得到新的 knowledge posterior。

接下來我們我推導一件事情，就是我只知道 prior 是 multivariate gaussian，我也知道 likelihood 是 univariate gaussian，但我們知道他們相乘後出來的機率分佈仍是 multivariate？

接下來就來推導：
首先，在 likelihood 中的 a 先假設為常數，且每一個 x 所對應到的 variance 都一樣，稍後再說明原因。

因此，
${\frac{-a}{2}(||Xw-y||^2\ +\ \frac{b}{a}w^Tw)\ \propto\ a||Xw-y||^2\ +\ bw^Tw \\
=\ a(Xw-y)^T(Xw-y)\ +\ bw^Tw \\
=\ a(w^TX^TXw\ -\ 2w^TX^Ty\ +\ y^Ty)\ +\ bw^Tw \\
=\ w^T(aX^TX+bI)w\ -2aw^TX^Ty\ +\ ay^Ty}$

對應到 ==multivariate gaussian 的 quadratic form==：
${(w-\mu)^T\Sigma(w-\mu) \\
=\ w^T\Sigma w\ -\ 2w^T\Sigma\mu\ +\ \mu^T\mu}$

在做對照時不用考慮常數項，因此可以得到：
${\Sigma\ =\ aX^TX+bI \\
aw^TX^Ty\ =\ w^T\Sigma\mu\ \\
\Rightarrow\ \mu\ =\ \Sigma^{-1}aX^Ty\ =\ a\Sigma^{-1}X^Ty}$

因此可以說明，更新後的 posterior 仍可以表示成一個 multivariate gaussian：
${P(w|D)\sim N(a\Sigma^{-1}X^Ty,\ aX^TX+bI)}$

那如果==中心不是的 0 且 variance 彼此不相等==時呢？
${P(w|D)\ =\ e^{\frac{-a}{2}(y_i-x_iw)^2+\frac{-1}{2}(w-m)^TS(w-m)} \\
\begin{split}
\Rightarrow\ max\ &=\ \frac{-a}{2}(||Xw-y||^2\ +\ \frac{1}{a}(w-m)^TS(w-m))\ \propto\ a||Xw-y||^2\ +\ (w-m)^TS(w-m) \\
&=\ a(Xw-y)^T(Xw-y)\ +\ (w-m)^TS(w-m) \\
&=\ a(w^TX^TXw\ -\ 2w^TX^Ty\ +\ y^Ty)\ +\ (w^TSw\ -\ 2w^TSm\ +\ m^TSm) \\
&=\ w^T(aX^TX+S)w\ -2w^T(aX^Ty+Sm)\ +\ ay^Ty\ +\ m^Sm
\end{split}}$

一樣比較 quadratic form：
${w^T\Sigma w\ -\ 2w^T\Sigma\mu\ +\ \mu^T\mu}$

${\Sigma\ =\ aX^TX+S \\
w^T(aX^Ty+Sm)\ =\ w^T\Sigma\mu\ \\
\Rightarrow\ \mu\ =\ \Sigma^{-1}(aX^Ty+Sm)}$

<P style='page-break-after:always'></P>

接下來我們來看，在更新的過程， prior、likelihood 與 posteior 是怎麼被更新的。


在上圖中，我們假設初始的 prior mean 為 0，variance 為 1，這時我們有第一筆 data，因此我們做了 piror 與 likelihood 的更新，可以看到，在第二張 piror 中，${w}$ 的分佈被更新的，而左圖則是在某一點可能的 ${(w_0,\ w_1)}$ 下，likelihood 的分佈，==中心顏色最深的部份，就是我們所想找的 fit line==。

> 可以回想前面說的，==對於每一個 x，我們都想找一個 gaussian distribution 去 model 它==。

> 這裡順便解釋 a，在 likelihood 中的 a，其實==不太影響最後找到的最佳解==，它只是你希望每個 x 的那個 distribution 有多寬，所以如果接近零，就會得到一條直線，但在訓練的過程，你也會希望==給予這條線以外的其他地方一些機率==，所以會給予一個定值 (沒理由也沒必要每個 x 的 variance 不一樣)。

所以可以發現，當迭代到最後，==${w}$ 縮到一點==，那就是我們想到的參數，而觀察右圖的 data space，可以發現，當迭代的次數增加後，不同可能的 ${w}$ (高機率區) 所得到的 fit line 越來越相近。

> 與 rLSE 比較不同的地方在於，即使初始點選得不好，還是會往正確的方向移動，而不是有一個假定的固定值。

## Fully Bayesian

在 MAP 中，我們可以找到最合適的參數 w，但實際上，我們想知道的是從 x 到 y 的機率分佈，我們並不僅僅只是想知道在特定 w 下的 y 是長怎樣，也想知道其他的 w 下，y 的分佈是如何，也就是想要得到一個 x 到 y 的機率分佈。

而這其實就是 marginalize，也就是 ${P(Y|\theta)\ =\ \int P(Y|w, \theta)\ P(w|\theta)}$。

這裡可以發現，${P(Y|w, \theta)}$ 是 multivariate Gaussian，而 ${P(w|\theta)}$ 則是 univariate Gaussian，在 marginalize 後，得到的 Gaussian 仍為 multivarite （未證明）。

推導：
${\begin{split}令\ &P(Y|w, \theta)\ =\ N(Xw, a) \\
&P(w|\theta)\ =\ N(\mu,\Sigma)\end{split}\\
\begin{split}則\ P(Y|\theta)\ &=\ \int N(Xw, a)\ N(\mu,\Sigma)dw\ =\ \int e^{\frac{-a^{-1}}{2}(Xw-Y)^T(Xw-Y)}e^{\frac{-1}{2}(w-\mu)^T\Sigma^{-1}(w-\mu)} \\
&=\ \int e^{\frac{-1}{2}(a^{-1}(w^TX^TXw-2w^TX^TY+Y^TY)+w^T\Sigma^{-1} w-2w^T\Sigma^{-1}\mu+\mu^T\mu)} \\
&=\int e^{\frac{-1}{2}(w^T(a^{-1}X^TX+\Sigma^{-1})w-2w^T(a^{-1}X^TY+\Sigma^{-1}\mu)+a^{-1}Y^TY+\mu^T\mu)}\end{split}}$

和前面一樣，我們對照 quadratic form：
${w^TC^{-1} w\ -\ 2w^TC^{-1}m\ +\ m^TC^{-1}m}$

${\begin{split}\Rightarrow\ &C^{-1}\ =\ a^{-1}X^TX+\Sigma^{-1} \\
&m\ =\ C(a^{-1}X^TY+\Sigma^{-1}\mu)
\end{split}}$

所以指數項可以變成：
${\begin{split}
&w^TC^{-1} w\ -\ 2w^TC^{-1}m\ +\ m^TC^{-1}m\ -\ m^TC^{-1}m\ +\ a^{-1}Y^TY\ +\ \mu^T\mu \\
\Rightarrow\ &P(Y|\theta)\ =\ \int e^{w^TC^{-1} w\ -\ 2w^TC^{-1}m\ +\ m^TC^{-1}m}e^{-m^TC^{-1}m\ +\ a^{-1}Y^TY\ +\ \mu^T\mu}dw \\
&=\ e^{-m^TC^{-1}m\ +\ a^{-1}Y^TY\ +\ \mu^T\mu}\int e^{w^TC^{-1} w\ -\ 2w^TC^{-1}m\ +\ m^TC^{-1}m}dw\end{split}}$

又我們知道，Gaussian 的積分和為 1
${\Rightarrow\ P(Y|\theta)\ =\ e^{-m^TC^{-1}m\ +\ a^{-1}Y^TY\ +\ \mu^T\mu}}$

指數項帶入 ${m}$ 可得：
${\begin{split}&-a^{-2}Y^TXCX^TY\ -\ 2a^{-1}Y^TXC\Sigma^{-1}\mu\ +\ a^{-1}Y^TY\ +\ \mu^T\mu \\
=\ &Y^T(a^{-1}\ -\ a^{-2}XCX^T)Y\ -\ 2a^{-1}Y^TXC\Sigma^{-1}\mu\ +\ 常數\end{split}}$

${令\ \lambda\ =\ a^{-1}\ -\ a^{-2}XCX^T \\
\Rightarrow\ Y^T\lambda Y\ -\ 2a^{-1}Y^TXC\Sigma^{-1}\mu\ +\ 常數}$

再對照一次 quadratic form：
${Y^TC^{'-1}Y\ -\ 2Y^TC^{'-1}m^{'}\ +\ m^{'T}C^{'-1}m^{'}}$

${\begin{split}\Rightarrow\ &C^{'-1}\ =\ \lambda \\
&m^{'}\ =\ C^{'}(a^{-1}XC\Sigma^{-1}\mu)\ =\ \lambda^{-1}(a^{-1}XC\Sigma^{-1}\mu)
\end{split}}$

接著我們要把 ${C\ 與\ \lambda}$ 帶入，首先由 Sherman–Morrison formula 得：
若 ${C^{-1}\ =\ a^{-1}X^TX+\Sigma^{-1}}$
則 ${C\ =\ \Sigma\ -\ \frac{\Sigma a^{-1}X^TX\Sigma}{1+a^{-1}X\Sigma X^T}}$

為了方便計算，令 ${\alpha\ =\ X\Sigma X^T}$：
${C\ =\ \Sigma\ -\ \frac{\Sigma a^{-1}X^T\alpha X^{-T}}{1+a^{-1}\alpha}}$

因此，
${\begin{split}C^{'-1}\ &=\ \lambda\ =\ a^{-1}\ -\ a^{-2}XCX^T \\
&=\ a^{-1}\ -\ a^{-2}X(\Sigma\ -\ \frac{\Sigma a^{-1}X^T\alpha X^{-T}}{1+a^{-1}\alpha})X^T \\
&=\ a^{-1}\ -\ a^{-2}(\alpha\ -\ X(\frac{\Sigma a^{-1}X^T\alpha}{1+a^{-1}\alpha})) \\
&=\ a^{-1}\ -\ a^{-2}(\alpha\ -\ \frac{a^{-1}\alpha^2}{1+a^{-1}\alpha}) \\
&=\ a^{-1}\ -\ a^{-2}(\frac{\alpha}{1+a^{-1}\alpha}) \\
&=\ \frac{a^{-1}}{1+a^{-1}\alpha}\end{split} \\
\Rightarrow\ C{'}\ =\ \frac{1+a^{-1}\alpha}{a^{-1}}\ =\ a\ +\ \alpha\ =\ a\ +\ X\Sigma X^T}$

${\begin{split}而\ m^{'}\ &=\ \lambda^{-1}(a^{-1}XC\Sigma^{-1}\mu)\ =\ a^{-1}(\lambda^{-1}XC\Sigma^{-1}\mu)\ =\ a^{-1}(\lambda^{-1}XC\Sigma^{-1}\mu)^T \\
&=\ a^{-1}(\mu^T\Sigma^{-1}CX^T\lambda^{-1})\ =\ a^{-1}\mu^T\lambda^{-1}\Sigma^{-1}CX^T \\
&=\ a^{-1}\mu^T\lambda^{-1}\Sigma^{-1}(\Sigma\ -\ \frac{\Sigma a^{-1}X^T\alpha X^{-T}}{1+a^{-1}\alpha})X^T \\
&=\ a^{-1}\mu^T\lambda^{-1}(X^T\ -\ \frac{a^{-1}X^T\alpha}{1+a^{-1}\alpha})\ =\ a^{-1}\mu^T\lambda^{-1}X^T(\frac{1}{1+a^{-1}\alpha}) \\
&=\ \mu^T\lambda^{-1}X^T(\frac{a^{-1}}{1+a^{-1}\alpha})\ =\ \mu^T\lambda^{-1}X^T\lambda\ =\ \mu^TX^T\ =\ X\mu\end{split}}$

所以可以得到，${P(Y|\theta)\ =\ N(X\mu, a+X\Sigma X^T)}$

由上面的式子可以發現，變異數由兩項相得所得，當 MAP 的迭代次數增加時，其 ${\Sigma}$ 會趨近於零，也就是我們對所得到的 ${w}$ 有一定的信心，變異數就只剩下 ${a}$ 了。