從貝氏的世界看 curve fitting: Maximum Likelihood (ML)、Maximum a Posteriori (MAP)

今假設我們有一組資料數據，$\mathbf{X} = \{\mathbf{x}_1, \dots, \mathbf{x}_N\}^\top$ 是這組資料的輸入值，而 $\mathbf{t} = \{t_1, \dots, t_N\}$ 是這組資料的輸出值，也是我們想要估計的對像。其中 $\mathbf{x}_n$ 為一個 column vector，此向量裡的每個元素都代表不同的意義。舉例來說，$\mathbf{t}$ 可以是台積電(2330)上市以來的每日收盤價，而 $\mathbf{x}_n$ 可能對應的是第 $n$ 日的 EPS、營收、KD 值等等的數據，我們想要利用既有的 $\mathbf{X}$ 以及 $\mathbf{t}$ 資料數據建構一個可以解釋輸入與輸出關係的模型，一旦有了此模型，我們就可以透過未來的台積電資訊 $\mathbf{x}_{N+1}$ 估計出未來台積電可能的價位。

貝氏的世界裡將輸入值和輸出值的關係看作一個 deterministic 函式 $y(\mathbf{x}_n, \mathbf{w})$ 搭配上一個 IID 的 Gaussian noise $\epsilon = \mathcal{N}(\epsilon | 0, \beta^{-1})$:

$$t_n = y(\mathbf{x}_n, \mathbf{w}) + \epsilon$$

其中 $\mathbf{w}$ 是函式 $y$ 中的未知參數。延續剛才股票模型的例子，假設我們預設的輸入值 $\mathbf{x}_n$ 跟股價 $t_n$ 真的有某種函式關係，現實生活中也會因為各種不可控制的變數影響最後的股價，例如投資者過度樂觀或是過度恐懼等等，都會導致股價沒有照著預期的模型走，而這些未知的、無法控制的因素我們都用 $\epsilon$ 這個 Gaussian 隨機變數來表示，因為根據中央極限定理 (Central Limit Theorem)，大量相互獨立的隨機變數之和趨近於常態分布，因此 Gaussian noise 是最適合作為未知因素的隨機變數。

接下來我們可以得到輸出值 $t_n$ 的機率分布：

$$p(t_n|x_n, \mathbf{w}) = p(t_n|x_n) = \mathcal{N}(t_n | y(\mathbf{x}_n, \mathbf{w}), \beta^{-1})$$

注意到我們在 maximum likelihood 裡會將 $\mathbf{w}$ 視為一個固定量，而非 random variable，因此可以得到 $p(t_n|x_n, \mathbf{w})$ 和 $p(t_n|x_n)$ 的分布其實相同，但是為了後續 maximum a posteriori 的推導，我們還是將 $\mathbf{w}$ 保留在條件機率裡面。另外，我們在這裡做出了一個假設，就是所有資料 $\mathbf{X} = \{\mathbf{x}_1, \dots, \mathbf{x}_N\}^\top$ 和 $\mathbf{t} = \{t_1, \dots, t_N\}$ 皆是獨立取樣出來的資料，並不會互相影響。因此他們的聯合分布可以寫為：

$$p(\mathbf{t}|\mathbf{X}, \mathbf{w}) = \prod_{n = 1}^N \mathcal{N}(t_n | y(\mathbf{x}_n, \mathbf{w}), \beta^{-1})$$

此時我們可以利用 maximum likelihood 的方法估計未知的模型參數 $\mathbf{w}$ 和 $\beta$，即是尋找某個 $\mathbf{w}_{\text{ML}}$ 以及 $\beta_{\text{ML}}$ 使得所有資料的 joint maximum likelihood 有最大值：

$$\mathbf{w}_{\text{ML}} = \arg \max_{\mathbf{w}} \prod_{n = 1}^N \mathcal{N}(t_n | y(\mathbf{x}_n, \mathbf{w}), \beta^{-1})$$ $$\beta_{\text{ML}} = \arg \max_{\beta} \prod_{n = 1}^N \mathcal{N}(t_n | y(\mathbf{x}_n, \mathbf{w}), \beta^{-1})$$

由於上述的數學式為一個凸函式(convex function)，因此我們可以直接利用微分等於 $0$ 的特性求解。為了使計算微分更容易，我們將聯合分布 $p(\mathbf{t}|\mathbf{X})$ 取 $\log$，再做運算。注意 $\log$ 是一個單調遞增(monotonically increasing)函式，並不會影響極值發生的位置，又可以使推導的過程變得更加容易：

$$\begin{aligned} &p(\mathbf{t}|\mathbf{X}, \mathbf{w}) = \frac{1}{\sqrt{(2\pi)^k |\beta^{-1} \mathbf{I}|}} \exp \bigg[-\frac{1}{2} \big(\mathbf{t} – y(\mathbf{X}, \mathbf{w}) \big)^\top (\beta^{-1} \mathbf{I})^{-1} \big(\mathbf{t} – y(\mathbf{X}, \mathbf{w}) \big) \bigg] \\ \Rightarrow &\ln p(\mathbf{t}|\mathbf{X}, \mathbf{w}) = -\frac{\beta}{2} \sum_{n = 1}^N \{y(\mathbf{x}_n, \mathbf{w}) – t\}^2 + \frac{N}{2} \ln \beta – \frac{N}{2} \ln (2\pi) \end{aligned}$$

接下來我們就可以計算 $\ln p(\mathbf{t}|\mathbf{X})$ 對 $\mathbf{w}$ 和 $\beta$ 的微分，求出 $\mathbf{w}_{\text{ML}}$ 以及 $\beta_{\text{ML}}$，為避免篇幅過長，將這部分的詳細計算將留到Machine Learning 入門：多項式擬合 (Polynomial Curve Fitting) 與 regularization 推導。值得注意的是，找到某個 $\mathbf{w}_{\text{ML}}$ 使得 log-likelihood 有最大值，其實等同於找到某個 $\mathbf{w}_{\text{ML}}$ 使得 $\frac{1}{2} \sum_{n = 1}^N \{y(\mathbf{x}_n, \mathbf{w}) – t\}^2$ 有最小值。可以直接忽略 $\beta$ 是因為不論其值為何都不會影響最後 $\mathbf{w}_{\text{ML}}$ 的結果。從這裡我們得到一個有趣的結論，利用 Gaussian noise 的機率模型建構出來的 maximum likelihood，其極值發生的位置 $\mathbf{w}_{\text{ML}}$ 剛好和 minimize mean square error 極值發生的位置相同！

接下來我們將 Maximum likelihood 延伸至 Maximum a Posteriori (MAP)，如同Maximum Likelihood & Maximum a Posteriori－基礎估計模型的詳細介紹裡面提到的，在現實中我們可能常常因為收集到的 dataset 有某些程度的 bias，導致訓練出來的模型效果不如預期。為了避免模型過多依賴手頭上的 dataset，我們會多加一個 prior distribution，藉此將我們對整個問題的先驗知識導入在模型當中。舉例來說，假設我們收集到的台積電股市資料集非常偏差，股價完全等於 30 倍本益比 (本益比 = 股價 / EPS)，那我們估計出來的 $\mathbf{w}$ 權重可能完全偏向 $\mathbf{x}_n$ 裡面的 EPS，其他輸入如 KD 值、營收等等就完全不列入考慮了，深具經驗的投資人如果看到這個模型肯定會提出疑惑。為了避免 dataset 裡面的偏差過多影響模型的結果，我們可以基於對台積電的了解，針對 $\mathbf{w}$ 提出某種機率分布，例如我們可能認為營收比 EPS 來的重要，就可以假設 $\mathbf{w}$ 是一個 mean 不為 $0$ 的高斯分布，其對應到營收權重的 mean 可能稍微高一點，對應到 EPS 的 mean 可以稍微低一點。

我們想要利用 posterior 來估計未知的 $\mathbf{w}$，就必須推導出 $p(\mathbf{w} | \mathbf{X}, \mathbf{t})$： $$ p(\mathbf{w} | \mathbf{X}, \mathbf{t}) = \frac{p(\mathbf{w}, \mathbf{X}, \mathbf{t})}{p(\mathbf{X}, \mathbf{t})} = \frac{p(\mathbf{t} | \mathbf{X}, \mathbf{w}) p(\mathbf{X})p(\mathbf{w})}{p(\mathbf{X}, \mathbf{t})} \propto p(\mathbf{t} | \mathbf{X}, \mathbf{w})p(\mathbf{w})$$

這裡我們做出一個假設，即 $\mathbf{X}$ 和 $\mathbf{w}$ 為獨立分布。因為輸入值如何生成、改變，其實都跟我們建構的模型沒有關係，台積電的 EPS、本益比，也不會被我們的權重 $\mathbf{w}$ 是多少而影響。或是我們可以用最直觀的方式看出兩個隨機向量互相獨立：我們已經將 $\mathbf{w}$ 假設成一個 zero mean 且 variace 為 $\alpha$ 的 Gaussian 了，那其分布和 $\mathbf{X}$ 獨立也就顯而易見了。另外，因為 $\mathbf{X}$ 以及 $\mathbf{t}$ 對於我們來說是已知的固定量，$p(\mathbf{X}, \mathbf{t})$ 的值是多少也是已知的 (縱使我們沒有對該分布做出任何假設，但我們還是把它當作已知)，因此當我們在估計 $\mathbf{w}$ 時，可以忽略 $p(\mathbf{X}, \mathbf{t})$ 這項，直接利用 $p(\mathbf{t} | \mathbf{X}, \mathbf{w})p(\mathbf{w})$ 來估計 $\mathbf{w}$。

我們遵循 ML 的方式，先取 log 再求最大值，可以推導出： $$\begin{aligned} &p(\mathbf{t}|\mathbf{X}, \mathbf{w}) p(\mathbf{w}) = \frac{1}{\sqrt{(2\pi)^k |\beta^{-1} \mathbf{I}|}} \exp \bigg[-\frac{1}{2} \big(\mathbf{t} – y(\mathbf{X}, \mathbf{w}) \big)^\top (\beta^{-1} \mathbf{I})^{-1} \big(\mathbf{t} – y(\mathbf{X}, \mathbf{w}) \big) \bigg] \frac{1}{\sqrt{(2\pi)^k |\alpha^{-1} \mathbf{I}|}} \exp \bigg[-\frac{1}{2} \mathbf{w}^\top (\alpha^{-1} \mathbf{I})^{-1} \mathbf{w} \bigg] \\ \Rightarrow &\ln \{p(\mathbf{t}|\mathbf{X}, \mathbf{w})p(\mathbf{w})\} = -\frac{\beta}{2} \sum_{n = 1}^N \{y(\mathbf{x}_n, \mathbf{w}) – t\}^2 + \frac{N}{2} \ln \beta – \frac{N}{2} \ln (2\pi) – \frac{\alpha}{2} \mathbf{w}^\top \mathbf{w} \end{aligned}$$