【翻譯】擬合與高斯分佈 [Curve fitting and the Gaussian distribution]

參考與前言

1. 英文原版 Original English Version：https://fabiandablander.com/r/Curve-Fitting-Gaussian.html
2. 原文中有超多參考，原文參考就不一一複製過來了哈

注意：此文主要是 英文原版的自制中文翻譯 (已徵得原作者同意)，而且可能進行了一定的壓縮/新增操作，請學有餘力者跳轉至原版連結

Note: This post was originally published on  https://fabiandablander.com/

沒錯，這裡是我看大佬的傑作，一開始提出的問題，其實類似於回顧吧 回顧為什麼之旅：為什麼在擬閤中是高斯分佈？高斯過程又在幹什麼？而這篇主要是針對前者的回覆，擬合和高斯分佈 (後者關於高斯過程中的下次再搜一搜) 然後這篇英文文章真的超讚！作者還有其他系列的也超讚！主要是看的過程還可以學歷史?hhhh 由淺入深，看看我有沒有精力再把第二篇關於高斯的翻完 (原文地址：https://fabiandablander.com/statistics/Two-Properties.html) 讀完後覺得也很贊

省流版：中心極限定理

擬合的過程

---

一條線在平面座標系的表示主要是斜率(\(b_1\))和截距(\(b_0\))（也就是x=0時，y的那個點）

\[y=b_0+b_1x \]

那麼假設我們擁有很多個點後去得到 \(b_1\) 和 \(b_0\) 。先從簡單的來吧

1. 假設只有一個點的時候，我們的擬合結果可能像第一幅，也就是無數條線，無法確定 \(b_1\) 和 \(b_0\)​
2. 當我們有兩個點時 我們可以唯一確定下來這條直線是誰
3. 當我們有三個點時，我們無法找到一條直線把三個點都穿過，但是可以想想，是否畫出來的線中 有較好的線(擬合誤差較小)

那麼接下來，讓我們繼續思考，第三幅圖我們應該怎樣去擬合？或許是再加一個引數 \(b_2\)​ ？那 \(b_2\)​ 乘什麼呢？還是簡單起見，我們就直接按 \(x\)​ 的指數上去把給他乘一個 \(x^2\)​ ，那現線上方程變成了這樣：

\[y=b_0+b_1x+b_2x^2 \]

現在是否發現了一些事情，在絕對的情況下，兩個引數 \((b_0,b_1)\)​ 用兩個點 \(P_1,P_2\)​ ，三個引數\((b_0, b_1, b_2)\)​ 用三個點 \(P_1,P_2, P_3\)​​

另外，如果我們同時可以這麼去看點 \(P'=(y,x_1,x_1^2)\)​ ，還是第三幅圖的點座標，我們可以用線性代數裡的矩陣形式去表示等式：

\[\begin{aligned} \mathbf{y} &=\mathbf{X} \mathbf{b} \\ \left(\begin{array}{l} 2 \\ 1 \\ 2 \end{array}\right) &=\left(\begin{array}{ccc} 1 & 1 & 1 \\ 1 & 3 & 9 \\ 1 & 2 & 4 \end{array}\right) \cdot\left(\begin{array}{l} b_{0} \\ b_{1} \\ b_{2} \end{array}\right) \end{aligned} \]

而這也就是上面的擬合結果，然後我們再思考一下：這樣的結果有什麼缺點？

1. 過擬合：對現在已知的資料點擬合完美，但是對未來的資料點可能就不OK了

2. 這個資料的意義是什麼：we haven’t really explained anything

   引用Fisher的話:

   The objective of statistical methods is the reduction of data. A quantity of data, which usually by its mere bulk is incapable of entering the mind, is to be replaced by relatively few quantities which shall adequately represent […] the relevant information contained in the original data. (Fisher, 1922, p. 311)

總結起來就是：沒錯 我們擬合了這條線，但是我們並沒有降低資料量，三個點還是需要對應的三個引數，那是否n個點就需要對應的n個引數呢？所以引出了回答上面的高亮部分：哪條直線更好？

"最好"的擬合

---

將超定(overdetermined)減小到確定(determined)情況並不可以。但是可以將超定轉到欠定。為了達到這一目的呢，我們需要做出合理的假設說明觀測點（資料點）是有對應的誤差的，例如這樣子：

\[\begin{aligned} &1=b_{0}+4 b_{1}+\epsilon_{1} \\ &4=b_{0}+2 b_{1}+\epsilon_{2} \\ &3=b_{0}+1 b_{1}+\epsilon_{3} \end{aligned} \]

而這裡的 \((\epsilon_{1}, \epsilon_{2}, \epsilon_{3})\)​​​​ 是unobserved quantities 未知誤差值。而這一步沒錯又引入了未知量，所以我們再次以少的資料（少的方程）去解多的引數量，所以我們的之前的超定就變成了欠定系統求解了

正如前面所說，只有一個點時，如果沒有其他限制，我們就有無數條直線說 我們擬合成功了。一樣的道理在這裡，我們需要新增限制，例如前面例子我需要經過這個點的線同時儘可能平行x軸，這就是一個限制。而對於現在這個問題 法國數學家Adrien-Marie Legendre 提出了著名的 最小二乘法

在上面那個例子中，通過最小二乘的誤差限制，只會出現一條線滿足最小的誤差。因此，對於所有的其他欠定系統我們都可以這樣：新增最小二乘誤差的限制，例如下圖：

同時最小二乘法的發展也是數理統計的分水嶺——Stephen Stigler 將其重要性比作數學中微積分的發展。那麼理論上我們已經知道可以這樣去找到最好的擬合。那數學上我們怎麼知道是具體哪一條，這裡提到了兩種，一種是由Legendre提出的優化；另外一種是幾何意義上的直觀insight

Least squares I: Optimization

首先我們的目標是找到最小的均方誤差。從簡單入手呢，就是直接求資料的平均值，然後再去減他們，類似於這樣，這個步驟也叫中心化 (centered the data)：

\[\begin{aligned} y'&=y-\frac{1}{n}\sum_{i=1}^n y_i \\ x'&=x-\frac{1}{n}\sum_{i=1}^n x_i \end{aligned} \]

這樣子的預估方式會使得 \(b_0=0\) ，所以就省略了去預估 \(b_0\)​ 是多少了 (未知量 bi--音效 減少一個)；原因是：中心化後資料平均值為0，那麼我們就把所有點相加：(懶得打公式了 我手寫了)

注意，下面的描述 \(x,y\) 都是已經中心化後的結果

1. 首先我們觀測到資料點 \(y_i\) 我們直線給出的預估是 \(x_ib_1\) ，所以誤差就是： \(\epsilon_{i}=y_{i}-x_{i} b_{1}\)​，我們把有的所有資料都加起來 那總誤差就是這樣的：

   \[\sum_{i=1}^{n} \epsilon_{i}^{2}=\sum_{i=1}^{n}\left(y_{i}-x_{i} b_{1}\right)^{2} \]

2. 現在我們就去找哪一個 \(b_1\) 可以得到最小的總誤差，這裡用 \(\hat b_1\) 表示

   \[\hat{b}_{1}=\underbrace{\operatorname{argmin}}_{b_{1}}\left(\sum_{i=1}^{n}\left(y_{i}-x_{i} b_{1}\right)^{2}\right) \]

3. 就像我們找曲線的最低點一樣，求導，導數等於0，值小於其周圍值的地方就是我們想要的點。那麼在這個問題中，我們就是這樣子做：

   \[\begin{aligned} 0 &=\frac{\partial}{\partial b_{1}}\left(\sum_{i=1}^{n}\left(y_{i}-x_{i} b_{1}\right)^{2}\right) \\ 0 &=\frac{\partial}{\partial b_{1}}\left(\sum_{i=1}^{n} y_{i}^{2}-2 \sum_{i=1}^{n} y_{i} x_{i} b_{1}+\sum_{i=1}^{n} x_{i}^{2} b_{1}^{2}\right) \\ 0 &=0-2 \sum_{i=1}^{n} y_{i} x_{i}+2 \sum_{i=1}^{n} x_{i}^{2} b_{1} \\ 2 \sum_{i=1}^{n} x_{i}^{2} b_{1} &=2 \sum_{i=1}^{n} y_{i} x_{i} \\ \hat{b}_{1} &=\frac{\sum_{i=1}^{n} y_{i} x_{i}}{\sum_{i=1}^{n} x_{i}^{2}} \end{aligned} \]

   其中 \(\sum_{i=1}^n y_i x_i\) 是\(x,y\) 之間的協方差 (scaled by \(n\)), and \(\sum_{i=1}^n x_i^2\) 是\(x\)自身的方差

Least squares II: Projection

從其他角度去看這個最好的擬合線，就是幾何了，首先，直接把點投影到擬合線上，計算這個y軸上的差距作為誤差 這個誤差是垂直於x軸的，就像下面這幅圖一樣

\[\begin{aligned} \left(\begin{array}{c} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{n} \end{array}\right) &=\left(\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array}\right)-\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ z_{n} \end{array}\right) b_{1} \\ & \epsilon=\mathbf{y}-\mathbf{x} b_{1} \end{aligned} \]

然後我們所希望的就是垂直於x軸的這個誤差很小，用數學來表示呢就是點乘加起來等於0，簡單版就是 \(\mathbf{x}^T \mathbf{\epsilon} = 0\)​​

\[\left(\begin{array}{llll} x_{1} & x_{2} & \ldots & x_{n} \end{array}\right)\left(\begin{array}{c} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{n} \end{array}\right)=0 \]

然後就又是最小二乘的方程式了：

\[\begin{aligned} \mathbf{x}^{T} \epsilon &=0 \\ \mathbf{x}^{T}\left(\mathbf{y}-\mathbf{x} b_{1}\right) &=0 \\ \mathbf{x}^{T} \mathbf{x} b_{1} &=\mathbf{x}^{T} \mathbf{y} \\ b_{1} &=\frac{\mathbf{x}^{T} \mathbf{y}}{\mathbf{x}^{T} \mathbf{x}} \\ b_{1} &=\frac{\sum_{i=1}^{n} x_{i} y_{i}}{\sum_{i=1}^{n} x_{i}^{2}} \end{aligned} \]

咦 這是不是長的很像我們在方法一 優化裡得到的結果，沒錯就是一模一樣的。

最後，可以注意到我們一直沒有對 \(b_0\) 進行估計，知道最小二乘結束後 \(b_0\)​ 還是未知的，所以最後我們處理的方式是 \(\bar y = b_0\)​ 也就是 \(y\)​ 的平均值。也正是因為這一事實，高斯曾經將高斯分佈證明是誤差的分佈。

“多好才是最好” 高斯, 拉普拉斯

---

雖然呢，最小二乘給出了怎樣去擬合出“最好”的曲線來縮小誤差，但是他並沒說明

1. 誤差 \(\epsilon\)​​ 的隨機性
2. “多好才是最好”​ How good is best?

1809年是高斯將最小二乘放在了概率的角度去考慮。他做出了這樣一個假設：每一個誤差項 \(\epsilon_i\) 都是服從某個分佈 \(\phi\) 的
利用這分佈，當 \(\epsilon_i\) 很小時，\(\epsilon_i\)​​ 分佈整體概率很大，而這個情況也就是觀測和預測值十分接近的時候。然後我們再進一步假設所有誤差分佈都是 獨立同分布 ，他想找到能夠最大化概率 即最小化整體誤差

\[\Omega=\phi\left(\epsilon_{1}\right) \cdot \phi\left(\epsilon_{2}\right) \cdot \ldots \cdot \phi\left(\epsilon_{n}\right)=\prod_{i=1}^{n} \phi\left(\epsilon_{i}\right) \]

所以呢，現在我們就要去找到這樣的一個分佈 $\phi $。接著高斯大哥就發現我可以對 $\phi $ 做出一些說明，比如他一定是對稱的而且在等於誤差值=0的時候概率最大。接著他做出假設：平均值處應該就是能夠最好概括總結 \(n\) 個觀測值的地方 \((y_1,\dots,y_n)\) 所以，他假設最大化整個 \(\Omega\) 應該就和最小化整個均方誤差的和是一個性質

• [x] 這裡有一點點不太能理解，He then assumed that the mean should be the best value for summarizing \(n\)​ measurements

  作者在註釋裡寫了 the mean predicts the data best

• [ ] 這裡還有一個 should lead to the same solution as minimizing the sum of squared errors when we have one unknown
  為啥是當我們有一個未知時？是指 \(b_0\) 嘛？

然後繞完這個圈後，高斯就證明最小二乘的誤差分佈應該是這樣子的：

\[\phi\left(\epsilon_{i}\right)=\frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \left(-\frac{1}{2 \sigma^{2}} \epsilon_{i}^{2}\right) \]

其中 \(\sigma^{2}\)​ 是方差 (其實一開始高斯提出來的時候還不長這樣 最終為人所知的其實是Fisher修改了後方差後的版本) 這也就是為人所知的高斯分佈，雖然呢 他一開始是Stigler's law of eponomy，同時Karl Pearson也叫其正態分佈，雖然他對此有絲絲後悔：

“Many years ago I called the Laplace–Gaussian curve the normal curve, which name, while it avoids an international question of priority, has the disadvantage of leading people to believe that all other distributions of frequency are in one sense or another ‘abnormal’.” (Pearson, 1920, p. 25)

分佈圖如下：（雖然大家都知道了）

使用高斯分佈後，我們的最大化問題就變成了這樣子：

\[\Omega=\prod_{i=1} \phi\left(\epsilon_{i}\right)=\prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \left(-\frac{1}{2 \sigma^{2}} \epsilon_{i}^{2}\right)=\left(\frac{1}{\sqrt{2 \pi \sigma^{2}}}\right)^{n} \exp \left(-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{n} \epsilon_{i}^{2}\right) \]

注意：能使 \(\Omega\) 最大的值並不會隨其他的東西改變，比如添刪常數值，或者對其取對數log；

1. \(-\sum_{i=1}^{n} \epsilon_{i}^{2}\) 這是我們想要最大化的東西

2. 因為上面式子有負號，所以去掉負號的話是： \(\sum_{i=1}^{n} \epsilon_{i}^{2}\) ，同時最大化最小化這個

3. 最後再看一遍這個 \(\sum_{i=1}^{n} \epsilon_{i}^{2}\)​ 也正是最小化均方誤差！

   我們並沒有因為分佈遠離我們的初衷，到頭來分佈概率還在朝著同一個目標

Pierre Simone de Laplace 沒錯就是拉普拉斯大哥在1810年注意到了高斯分佈，然後給高斯誤差曲線再一次的證明：如果我們將誤差視為許多（微小）的擾動影響集合，這個集合將按照 中心極限定理正態分佈

中心極限定理

我覺得這裡直接源自百度百科可能更快：

中心極限定理，是指概率論中討論隨機變數序列部分和分佈漸近於正態分佈的一類定理。這組定理是數理統計學和誤差分析的理論基礎，指出了大量隨機變數近似服從正態分佈的條件。

前文中我們提到了散落的點是獨立同分布，在這個基礎上我們加的是 隨機性 ，也就是說有一串獨立同分布的隨機變數序列，如果此序列是有限的方差，那當 \(n\)​ 越大，此序列的平均值也就越大，而他也越接近正態分佈；當 \(n \rightarrow \infin\)​ 平均值實際收斂到了正態分佈

拉普拉斯意識到，如果將最小二乘問題中的誤差視為小影響的總和（即均值），那麼它們將呈正態分佈。 這為最小二乘解提供了an elegant justification

舉個python程式的例子吧（原文是R語言的哈）

首先我們假設由 \(m=500\) 個獨立同分布(程式中分佈為正態的) 生成的點，然後計算一個他們一起的誤差 \(\epsilon_i\)。假設我們有 \(n=200\) 組資料觀測，那麼就有200個獨立的誤差值，然後你可以可以發現這個誤差接近與高斯

import numpy as np
import matplotlib.pyplot as plt

n_error = 200
influence_one_error = 500

errors = list()

for i in range(n_error):
    errors.append(np.mean(np.random.uniform(-10, 10,influence_one_error)))

num_bins = 30

fig, ax = plt.subplots()

# the histogram of the data
n, bins, patches = ax.hist(errors, num_bins, density=True) # 誤差們直接的直方圖

# add a 'best fit' line
mu = 0
sigma = np.sqrt(20**2 / 12 ) / np.sqrt(influence_one_error) #隨意的 僅用來做表示
x = np.linspace(-1,1,num=1000)  # x軸的取值範圍
y = 1 / (sigma * pow(2 * np.pi, 0.5)) * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2))  # 概率密度函式公式

ax.plot(x, y, 'r--')
ax.set_ylabel('Probability density')
ax.set_title('Central Limit Theorem')

# Tweak spacing to prevent clipping of ylabel
fig.tight_layout()
plt.show()

生成的圖：

可能你還沒意識到這個有多棒，讓我幫你再來一次。首先我們是隨意的生成了500個正態分佈(按定義來講 其他分佈也OK 只要是他們一同屬於某一分佈) 的資料，然後算平均值，計算他們自己內部的誤差是多少；而我們有200組這樣的資料，那麼我們就收集到了每組資料中自己的誤差，然後我們發現誤差間是屬於正態分佈的。

線性迴歸

最後我們再回到線性迴歸這個問題，也就是之前提到的擬合。

高斯分佈的一個特點是任何正態分佈隨機變數的線性組合本身都是正態分佈的。 我們可以將線性迴歸問題寫成矩陣形式，這表明 y 是 x 的加權線性組合。 具體來說，如果我們有 n 個資料點，我們就有了一個由 n 個方程組成的系統，我們可以更簡潔地用矩陣表示法來寫

\[\begin{aligned} \left(\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array}\right) &=\left(\begin{array}{cc} 1 & x_{1} \\ 1 & x_{2} \\ \vdots & \vdots \\ 1 & x_{n} \end{array}\right) \cdot\left(\begin{array}{c} b_{0} \\ b_{1} \end{array}\right)+\left(\begin{array}{c} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{n} \end{array}\right) \\ \mathbf{y} &=\mathbf{X b}+\epsilon \end{aligned} \]

因為這種線性的關係，使得對誤差假設是正態分佈，同時也為 \(\mathbf{y}\) 本身帶來了這一假設條件，也就是服從正態分佈：

\[y_{i} \mid \mathbf{x}_{i} \sim \mathcal{N}\left(\mathbf{x}_{i}^{T} \mathbf{b}, \sigma^{2}\right) \]

換句話說，每個點的 \(y_i\) 概率密度可以由這個公式給出：

\[\frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \left(-\frac{1}{2 \sigma^{2}}\left(y_{i}-\mathbf{x}_{i}^{T} \mathrm{~b}\right)^{2}\right) \]

也就是下面這幅圖：

直觀上來看的話就是：誤差的方差越小，擬合的越好

後話 (相關性與迴歸係數)

---

相關性並不意味著因果關係。 有些人認為線性迴歸是一種因果模型，但是不是這樣的。

為了證明這一點，我們可以將簡單線性迴歸中的迴歸係數與相關性聯絡起來——它們僅在標準化上有所不同。

假設以均值為中心的資料，Pearson correlation的定義是：

\[r_{x y}=\frac{\sum_{i=1}^{n} x_{i} y_{i}}{\sqrt{\sum_{i=1}^{n} x_{i}^{2} \sum_{i=1}^{n} y_{i}^{2}}} \]

注意這裡的相關性是對稱的，他並不取決於我們是x與y相關或是y與x相關。但是呢，在迴歸中，我們使用x去預測y是這樣的：

\[b_{x y}=\frac{\sum_{i=1}^{n} y_{i} x_{i}}{\sum_{i=1}^{n} x_{i}^{2}} \]

但是我們如果是用y去預測x值，那係數就會變成這樣：

\[b_{y x}=\frac{\sum_{i=1}^{n} y_{i} x_{i}}{\sum_{i=1}^{n} y_{i}^{2}} \neq b_{x y} \neq r_{x y} \]

然而，通過對資料進行標準化，即通過將變數除以各自的標準差，迴歸係數就成為了樣本相關性，即

\[\begin{aligned} \frac{\partial L}{\partial b_{x y}} &=\frac{\partial}{\partial b_{x y}} \sum_{i=1}^{n}\left(\frac{y_{i}}{\sqrt{\sum_{i=1}^{n} y_{i}^{2}}}-b_{x y} \frac{x_{i}}{\sqrt{\sum_{i=1}^{n} x_{i}^{2}}}\right)^{2} \\ &=\frac{\partial}{\partial b_{x y}}\left(\frac{\sum_{i=1}^{n} y_{i}^{2}}{\sqrt{\sum_{i=1}^{n} y_{i}^{2}}}-2 b_{x y} \frac{\sum_{i=1}^{n} y_{i} x_{i}}{\sqrt{\sum_{i=1}^{n} y_{i}^{2} \sum_{i=1}^{n} x_{i}^{2}}}+b_{x y}^{2} \frac{\sum_{i=1}^{n} x_{i}^{2}}{\sum_{i=1}^{n} x_{i}^{2}}\right) \\ &=0-2 \frac{\sum_{i=1}^{n} y_{i} x_{i}}{\sqrt{\sum_{i=1}^{n} y_{i}^{2} \sum_{i=1}^{n} x_{i}^{2}}}+2 b_{x y} \\ 2 b_{x y} &=2 \frac{\sum_{i=1}^{n} y_{i} x_{i}}{\sqrt{\sum_{i=1}^{n} y_{i}^{2} \sum_{i=1}^{n} x_{i}^{2}}} \\ b_{x y} &=\frac{\sum_{i=1}^{n} y_{i} x_{i}}{\sqrt{\sum_{i=1}^{n} y_{i}^{2} \sum_{i=1}^{n} x_{i}^{2}}} \end{aligned} \]

可以看到最後他又等於了 \(r_{xy}\) ，同時 這個標準化的迴歸係數也可以通過乘以原始迴歸係數來實現，即

\[b_s = b_{xy} \times \frac{\sqrt{\sum_{i=1}^n x_i^2}}{\sqrt{\sum_{i=1}^n y_i^2}} = \frac{\sum_{i=1}^n y_i x_i}{\sum_{i=1}^n x_i^2} \times \frac{\sqrt{\sum_{i=1}^n x_i^2}}{\sqrt{\sum_{i=1}^n y_i^2}} = \frac{\sum_{i=1}^n y_i x_i}{\sqrt{\sum_{i=1}^n y_i^2 \sum_{i=1}^n x_i^2}} \]

同樣的對 \(b_{yx}\) 也是一樣的

\[b_s = b_{yx} \times \frac{\sqrt{\sum_{i=1}^n y_i^2}}{\sqrt{\sum_{i=1}^n x_i^2}} = \frac{\sum_{i=1}^n y_i x_i}{\sum_{i=1}^n y_i^2} \times \frac{\sqrt{\sum_{i=1}^n y_i^2}}{\sqrt{\sum_{i=1}^n x_i^2}} = \frac{\sum_{i=1}^n y_i x_i}{\sqrt{\sum_{i=1}^n y_i^2 \sum_{i=1}^n x_i^2}} \enspace \]