Formal Formulation?of Statistical Learning

統(tǒng)計學(xué)習(xí)和傳統(tǒng)的統(tǒng)計問題不太一樣。傳統(tǒng)的數(shù)理統(tǒng)計想要通過數(shù)據(jù)推斷模型參數(shù)，統(tǒng)計學(xué)習(xí)想要推斷的是函數(shù)。采用的數(shù)據(jù)是函數(shù)的input和output：{(x_i,y_i)}（這屬于supervised learning；沒有y_i屬于unsupervised learning）。

"A problem of finding a desired dependence using a?limited?number of observations."

這個說法有點含糊其辭。我們需要formal formulation。

G是一個input的generator，可能是隨機的，也可以是人為設(shè)計的。S是一個supervisor（有輸出，所以是監(jiān)督學(xué)習(xí)），相當(dāng)于一個黑箱，或者說它就是statistical model。它根據(jù)某個條件概率F(y|x)輸出隨機的y。F(y|x)是我們不知道的。LM是一個learning machine。它的取值在{f_\alpha(x)}這個函數(shù)類里面，通過training data{(x_i,y_i)}decide出一個f_\alpha(x)，作為supervisor的一個近似（supervisor本身未必有一個函數(shù)形式，但是可以想辦法用函數(shù)去近似它，比如后驗眾數(shù)等等）。

有了LM，我們有兩個目的：prediction（比如預(yù)測x_0輸入的輸出）與inference（比如比較兩個因素的效果哪個強）。

既然是要推斷函數(shù)，就有兩種：參數(shù)（限定函數(shù)空間）或者非參數(shù)（無窮維的函數(shù)空間上干事）。本文只考慮兩種簡單的參數(shù)問題，即限定為線性函數(shù)（linear regression）和indicator function（classification）。

仍然可以用傳統(tǒng)數(shù)理統(tǒng)計的哲學(xué)去思考：

1. 要一致地考慮統(tǒng)計模型的測度族中的所有參數(shù)；

2. 用decision theory的角度去看。拿到數(shù)據(jù)之后，根據(jù)數(shù)據(jù)構(gòu)建（decide）一個statistic （estimator/learning machine）\hat{f}（它是一個隨機函數(shù)而不是隨機變量了，歸于更廣泛的隨機元的定義，所以這個statistic得打上引號），這個statistic應(yīng)當(dāng)與實際的f比較接近。度量接近程度，仍然可以用傳統(tǒng)的MSE。但是現(xiàn)在有一個問題，statistic是隨機函數(shù)，不是一個數(shù)，沒法直接算MSE。我們可以考慮risk functional

來作為decision與真實值的偏差。把這個risk最小化的\alpha就是learning的“標準答案”。但是我們不知道F(x,y)，所以沒法得到這個確切的標準答案（其實我們已經(jīng)知道，如果是MSE的話，最優(yōu)的“標準答案”就是后驗均值，等等；但這是一個理想化的東西，我們只能盡量接近它）。統(tǒng)計學(xué)習(xí)要做的是minimizing the risk functional?on the basis?of?empirical data。所以要在某些點上算這個隨機函數(shù)的取值。這個“某些點”就會有講究了。如果我們在training data set里面取，得到的是training MSE，它并不是我們感興趣的，因為因為我們本來就知道這些點上面的取值，沒有prediction的必要了。我們感興趣的是test MSE，即在新的test data上函數(shù)的取值。評判一個statistical learning method的優(yōu)劣，用的應(yīng)該是test MSE。對于classification問題，上面的討論用training error rate和test error rate代替即可。

跟傳統(tǒng)統(tǒng)計學(xué)相比，感覺統(tǒng)計學(xué)習(xí)的方法都是highly practical的，并不是很在意統(tǒng)計模型到底咋樣，不管怎樣都行，直接學(xué)習(xí)就完事兒了。

一些trade-off

數(shù)理統(tǒng)計中出現(xiàn)的bias-variance trade-off現(xiàn)在在統(tǒng)計學(xué)習(xí)中有很多表現(xiàn)：

這些trade-off中根本性的一個是bias-variance trade-off。我們熟知MSE的分解公式：

理想的情況下，Var和bias應(yīng)該同時很小，但是這一點是無法做到的。在傳統(tǒng)的數(shù)理統(tǒng)計中，我們找不到一個optimal的estimator，會轉(zhuǎn)而考慮一些compromise，譬如Pareto optimal的estimator，比如限定在unbiased estimator中的UMVUE。但是unbiased并不是什么非得滿足的事情，雖然它看起來的確挺好，但是有時犧牲unbiasedness會換取更小的MSE。經(jīng)典的例子就是Shinkage estimator?，F(xiàn)在在統(tǒng)計學(xué)習(xí)中類似，var和bias按下葫蘆浮起瓢：

當(dāng)模型的flexibility上升（比如說自由度變多，變成更高階的多項式模型），bias當(dāng)然更小，但是隨機函數(shù)\hat{f}對于數(shù)據(jù)更加敏感，所以var更大。總的trade-off的結(jié)果就是，MSE（test MSE）呈現(xiàn)U-shape。Flexibility過高稱為過擬合，過低稱為欠擬合，U-shape最低點處才是最好的擬合。

上面討論的是test MSE。對于training MSE，當(dāng)然是flexibility越高越小的，overfitting越厲害，training MSE越小。不過我們并不看training MSE，而關(guān)注test MSE。

另外一個同時存在的trade-off就是解釋性的問題。如果flexibility很低，比如線性模型，那就很容易解釋，比如Y受X_1和X_2哪個影響更大，看看參數(shù)就知道了。但是flexibility比較高時，上面這種解釋性就削弱了。

線性回歸

線性回歸是所有統(tǒng)計模型中用得最多的。它是有明確的統(tǒng)計模型假定的（即supervisor的假定）：

注意同方差的假設(shè)（homoscedastic）是必要的，沒有它就干不了活。一般還假設(shè)\epsilon~N(0,\sigma^2)，不過它有時候可以去掉（對于Gauss-Markov定理）。

Learning machine想干的就是estimate \beta_0,\beta_1,\sigma^2。這就回到了傳統(tǒng)的統(tǒng)計問題。得到估計值\hat{\beta_0},\hat{\beta_1}之后，我們怎么evaluate它們？Evaluation基于residuals:

根據(jù)residual畫出的residual plot是直觀觀察一個regression效果的好東西：

而總的擬合效果按照residual sum of squares來考察：

Gauss-Markov定理

對于前面的假設(shè)，best linear unbiased estimator?(BLUE)為最小二乘估計。

Linear指的是estimator為y_i的線性組合，best指的是（在無偏的時候）var最小。

這個定理不需要對\epsilon的具體分布有任何假設(shè)。

具體來說，\beta_0和\beta_1的最小二乘估計為

記法很簡單：\beta_1是關(guān)于y線性的，而且大概是y/x這個樣子。\beta_0使得擬合曲線剛好過數(shù)據(jù)點的重心。

另外誤差項的無偏估計為

n-2是很好理解的。方括號里面剝掉了\bar{x}和\bar{y}的兩個自由度。

另外值得一提的是，Gaussian假設(shè)下這個最小二乘估計正好也是MLE。

下面給出多元的情況。多元的線性回歸模型的statistical model為

Y仍然是一個數(shù)，但是X是一個行向量，\beta是一個列向量。這個表達式強調(diào)了，我們說的linear是\beta的linear，X內(nèi)部包含X^2,sinX之類的都是無所謂的(所以線性回歸用來擬合多項式、bsinx之類的形式完全是可以的，這點從名字來看有點誤導(dǎo)性。但是sinbx這種就不行了)。X中通常把第一項取為常數(shù)1，表示截距。

現(xiàn)在我們有了training data，data仍然可以寫成上面的方程：

但是此時Y是一個列向量，X是一個矩陣了。X和Y包含了所有的數(shù)據(jù)。MATLAB里做回歸也是按照這種標準格式來的。

最小二乘的結(jié)果為（此時只有一個統(tǒng)計量\beta，所以其實比一元情況更簡潔）：

它是Y的線性組合（符合BLUE）。記法：X^{-1}Y，插入兩個X^T項。X本身不一定是方陣，所以Y=X\beta不是一個well-posed的線性方程問題，沒法用X^{-1}Y求解。但是我們有X的Moore–Penrose廣義逆(X^TX)^{-1}X^T，這就是替代X^{-1}的最小二乘解。

注意，這個寫法默認了X^TX是可逆的，即X的各列是線性無關(guān)的。若不然，說明自變量中存在嚴格的線性關(guān)系（collinearity），此時回歸問題不是良定的（三維空間中，一條線上的數(shù)據(jù)顯然沒法擬合出一個平面），必須剔除掉這個線性關(guān)系。

另外誤差項的無偏估計為

統(tǒng)計量的分布

做假設(shè)檢驗和區(qū)間估計都需要知道上面的最小二乘統(tǒng)計量的分布。它們的分布總結(jié)在下面的定理中：（都是在正態(tài)假設(shè)下；直接寫多元的情況，更加簡潔）

其中，\beta的分布直接可以通過\beta的表達式推出來。\sigma^2的分布也是顯然的。二者的獨立性類比于Fisher的那個定理。（這個分布說白了就是Fisher那個定理的推廣）

有了這些分布就可以做假設(shè)檢驗和區(qū)間估計了。比如\beta_1的1-\alpha的置信區(qū)間為

預(yù)測

現(xiàn)在我們想要對于input x_*的輸出做出預(yù)測。不僅想要知道估計值，而且需要估計區(qū)間。

這個問題說白了就是想要估計參數(shù)\beta_0+\beta_1x_*。

為此，我們需要知道\beta_0+\beta_1x_*的分布：

從而1-\alpha的置信區(qū)間為

這樣一來統(tǒng)計學(xué)習(xí)的兩個目的：inference和prediction都完成了。

Assessing the extent to which the model fits the data

有兩種辦法：residual standard error與R^2。

Residual standard error就是\sigma的estimator:

它是一種lack of fit的刻畫，但是問題在于它取決于Y的scale。給出一個RSE為100，你并不好判斷它到底算大還是小。所以更常用的是R^2統(tǒng)計量。它定義為

這樣理解：TSS是Y的總的variability，其中線性回歸能夠解釋的部分是TSS-RSS，所以這個比例越大，說明能夠解釋的比例越大，擬合越好。

這個R^2有一些問題。它是作為

的一個estimator，但是這是一個biased的estimator。其中一部分原因在于商的期望本來就與期望的商有偏差；另一部分原因在于自由度。這個問題在multi-dimentional regression中會變得很嚴重。如果考慮自由度的修正，我們可以很自然地得到下面這種adjusted R^2:

這是一個less biased estimator。它一定小于原來的R^2，所以肯定小于1，不過可能會變成負數(shù)。在高維情況下，這種修正是必要的。

線性回歸可能出現(xiàn)的若干問題

問題1:非線性性?？梢詮膔esidual plot看出來：

右邊沒有明顯的pattern，所以線性關(guān)系差不多是對的。左邊有一個U型的pattern，所以需要考慮加入非線性項。

問題2:異方差性。它也可以通過residual plot看出來：

左邊這種就有明顯的異方差性了。我們知道Fisher-Behrens問題迄今還沒解決，所以我們不指望有對于異方差性的完美的解決方案。我們能做的只是妥協(xié)：用上凸函數(shù)，比如lnY，把數(shù)據(jù)做變換，使異方差性減小。僅此而已。

問題3:outlier。指的是某個output y離預(yù)測值非常遠。它也可以通過residual plot看出來：

對于outlier，如果我們知道它是誤差導(dǎo)致的，那直接扔掉就行了。不過如果是因為模型的缺陷導(dǎo)致的，就要仔細考慮了。

怎么判別outlier？上面只是看著像。一個定量的判別準則是，畫studentized residual plot，即除以RSE，在這樣的圖上，超過+-3的點概率為0.27%，可以認定大概就是outlier了。

Outlier對于估計值的影響其實不算大。比如下面這張圖中，有沒有去掉outlier，擬合的曲線基本沒有變化。但是outlier對于區(qū)間估計影響很大。你想，一個很大的residual把estimated sigma提得很高，估計的區(qū)間就會變的異常大。

問題4：high leverage points。相對于outliers，它指的是x離得比較遠的點。它對于擬合曲線的影響很大（為什么呢？后面會說）：

如何判斷一個點是否屬于high leverage point?我們可以用Leverage Statistic?；氐阶钚《朔ǎ?/p>

紅色的這個projection matrix H把Y的數(shù)據(jù)投影到擬合值（它顯然是一個冪等矩陣）。這個矩陣的意義在于，h_{ij}刻畫了i這個input對于第j個output的預(yù)測的影響大小。因為冪等性，

所以對角元h_{ii}就刻畫了第i個input對于所有output的預(yù)測的影響的總大小，因此它稱為第i個輸入的leverage。所謂的high leverage指的就是這個leverage statistic比較大。

下面具體研究這個H矩陣的性質(zhì)。首先，h_{ii}\geqslant h_{ii}^2，故h_{ii}\leqslant 1。其次，顯然有h_{ii}\geqslant 0。接著考慮H的跡：

也就是說，平均的leverage是(p+1)/n（這里用p+1是為了把p作為predictor的數(shù)目。在X矩陣里，列數(shù)比predictor數(shù)目要多一，對應(yīng)著截距）。作為一個經(jīng)驗法則，如果某一個input的leverage超過了平均leverage的兩倍，我們就認為它可能是一個high-leverage point。

對于簡單線性回歸，我們可以很容易算出此時的leverage statistic為

其含義是很直觀的。

這個地方需要多說一句：前面好像默認了“離其他input比較遠的點，對擬合的影響比較大”。這個說法是對的，因為leverage statistic的定義完全是按照對于擬合的影響來定義的，而算出來的結(jié)果則跟這個點離其他input的距離有關(guān)。所以這個聯(lián)系是存在的。

問題5:共線性。幾個predictor之間的強線性關(guān)系會導(dǎo)致X接近于不滿秩，這樣雖然還是能做擬合，但是擬合的誤差會很大。（可以想像，我們想要在三維空間中的若干數(shù)據(jù)點來找出一個平面，盡可能接近這些數(shù)據(jù)點。但是如果數(shù)據(jù)點位于一條直線上，就有無窮多個符合要求的平面存在，回歸問題變成ill-posed。不過另一種情況是，比如說，（正的）X和X^2，雖然完全是一一對應(yīng)，但是并不是線性關(guān)系，所以不會破壞回歸。這就是多項式回歸的道理。在三維空間中，一條曲線擬合出一個平面是可行的。）

如何檢測collinearity？一個辦法是看協(xié)方差矩陣。如果其中有一項很大，說明對應(yīng)的兩個predictor之間有較強的線性關(guān)系。但是它的局限在于，如果是三個predictor之間的線性關(guān)系，就沒法檢測出了。所以更好的辦法是計算VIF（variance inflation factor），即某個\beta_j的方差從單變量模型道full model的擴大倍數(shù)。一個經(jīng)驗規(guī)則是，VIF超過5-10就指示有collinearity。

VIF有一個很直觀的計算方法：

這里指的是用去掉X_j之外的所有predictor來擬合X_j得到的R^2。如果真的存在某種線性關(guān)系，這個R^2會接近于1，所以VIF會非常大。

解決collinearity的方法：1.去掉一個predictor；2.把collinear的變量結(jié)合成一個，譬如說做均值作為一個新的predictor。

MATLAB實現(xiàn)

以MATLAB自帶的關(guān)于汽車的數(shù)據(jù)為例。我們想要獲知汽車的mpg（miles per gallon）與horsepower以及weight的關(guān)系，假設(shè)這個關(guān)系的形式為mpg=a+b*horsepower+c*weight+d*horsepower*weight。

load carsmall

x1 = Weight;

x2 = Horsepower;?

y = MPG;

X = [ones(size(x1)) x1 x2 x1.*x2];%注意這里是按照標準的回歸格式寫的

[b,bint,r,rint,stats]=regress(y,X)

scatter3(x1,x2,y,'filled')

hold on

x1fit = min(x1):100:max(x1);

x2fit = min(x2):10:max(x2);

[X1FIT,X2FIT] = meshgrid(x1fit,x2fit);

YFIT = b(1) + b(2)*X1FIT + b(3)*X2FIT + b(4)*X1FIT.*X2FIT;

mesh(X1FIT,X2FIT,YFIT)

xlabel('Weight')

ylabel('Horsepower')

zlabel('MPG')

view(50,10)

hold off

輸出的結(jié)果為：

• 四個參數(shù)的估計值：60.71，-0.0102，-0.1882，0.0000385。

• 四個參數(shù)的95%置信區(qū)間[51.39,70.03]等等。

• error variance的估計為15.2363354389984。

• R^2為0.7742。

回歸的曲面如下圖所示。

分類

分類問題的formulation就是：有一個supervisor把輸入（按照一定隨機性）歸入不同的category里面。我們希望訓(xùn)練出LM來模擬這個supervisor。這個LM的結(jié)果應(yīng)該是input空間的一個indicator function，或者等價地說就是分類的decision boundary。

三種常用的classifier是邏輯回歸，線性判別分析，以及KNN。

邏輯回歸

邏輯回歸是有明確的統(tǒng)計模型假定的（即supervisor的假定）：對于輸入x，y為1的概率為e^{\beta_0+\beta_1x}/(1+e^{\beta_0+\beta_1x})，其余概率為0。然后我們的learning machine想要學(xué)習(xí)的就是\beta_0和\beta_1。學(xué)習(xí)方法用的是MLE，不過這個MLE沒有解析解，都是直接用計算機數(shù)值求解的。訓(xùn)練完成之后，我們不僅可以做分類，而且可以具體地知道對于輸入x，它有多大的概率屬于第一類，多大的概率屬于第二類（在模型正確的前提下）。

邏輯回歸適用于多元輸入以及大于兩個category的分類，但是后者通常不用邏輯回歸。

MATLAB中邏輯回歸包含在fitglm這一廣義線性回歸的函數(shù)里面，用predict函數(shù)做預(yù)測。

線性判別分析(LDA)

LDA的思想極其簡單，就是Bayes；Bayes中不知道的概率就用正態(tài)近似去估計（均值估計\mu，pooled方差估計\sigma^2）。這樣LDA就一句話講完了。多元的情況完全一樣。如圖是一個經(jīng)典的分類結(jié)果，分為三類，邊界都是線性的。

LDA對supervisor的假定就是這個正態(tài)性。

LDA也能算出P(Y|X)，只不過這是一個近似的“概率”，應(yīng)當(dāng)理解成一種indication。

除此之外有些需要多講的地方。以兩個類別的分類為例，默認的LDA是按照后驗概率是否大于0.5這個閾值來判別的。實際上為了特別的需求可以用其他閾值。找閾值可以用ROC曲線：

還需要多講的一點是LDA中一個關(guān)鍵的假設(shè)：不同X|Y分布的同方差性。這又是一個同方差性假設(shè)（統(tǒng)計中到處都是），當(dāng)然也是為了方便。在這個假設(shè)下，boundary就是線性的，這也是LDA這個名字的由來。

如果去掉這個假設(shè)，看起來并不會變得復(fù)雜多少。的確如此。如果對于每個Y各自估計方差，這個方法就叫做QDA（quadratic discriminant analysis），此時邊界是二次曲線。

LDA和QDA之間如何抉擇呢？QDA參數(shù)更多，flexibility更高，在bias-variance trade-off下可能更好也可能更差，看具體情況。比如同方差假設(shè)badly off，那么QDA當(dāng)然更好啦。

MATLAB上做LDA的函數(shù)是fitcdiscr，然后用predict函數(shù)做預(yù)測。

KNN

KNN是一個不假定特定統(tǒng)計模型的非參數(shù)方法（任意的supervisor），所以可以適用于非線性的奇奇怪怪的邊界。

KNN也能算出P(Y|X)，只不過這是一個近似的“概率”，應(yīng)當(dāng)理解成一種indication。

選取K是一個比較關(guān)鍵的問題。K比較小對應(yīng)于很大的flexibility，一般會過擬合，可以想見邊界會很復(fù)雜。K比較大則會欠擬合。比如下面這張圖：

隨著K的減小，training error會下降（K=1的時候就完全沒有training error），而test error會呈現(xiàn)U型。如下圖所示：

下面是MATLAB實現(xiàn)。fitcknn用于構(gòu)建分類器，而predict用于分類。

首先用一個supervisor生成數(shù)據(jù)：

X=[];

Y=cell(0);

figure(1)

hold on

for i=1:500

?? x=rand*2*pi;?

?? y=randn;

?? X=[X;x,y];

?? r=rand;

?? p=exp(5*(y-sin(x)))/(1+exp(5*(y-sin(x))));

?? if r<p

?? ? ? Y=[Y,'red'];

?? ? ? plot(x,y,'or')

?? else

?? ? ? Y=[Y,'blue'];

?? ? ? plot(x,y,'ob')

?? end

end

這個supervisor的分界線是y=sinx：

接下來構(gòu)建分類器：?Mdl = fitcknn(X,Y,'NumNeighbors',5);然后就可以用label = predict(Mdl,[1,1])來做分類。predict(Mdl,[1,2])給出'red'，是一個正確的分類。

為了更加直觀地顯示出boundary，可以用以下代碼：

k=1;

Mdl = fitcknn(X,Y,'NumNeighbors',k);

C=zeros(63,31);

for x=1:63

? ? for y=1:31

? ? ? ? RE=predict(Mdl,[x/63*2*pi,y/31*3-1.5]);

? ? ? ? if strcmp('blue',RE{1})==1

? ? ? ? ? ? C(x,y)=1;

? ? ? ? end

? ? end

end

pcolor(C')

然后我們改變k，看看邊界怎么樣。

k=1：過擬合非常嚴重

k=2:過擬合稍微好了一點

k=5:過擬合還是有

k=10:比較舒服了

k=20:看起來是最好的分類結(jié)果

k=50:欠擬合

k=100:嚴重的欠擬合

由此可以直觀地看出k對于KNN的影響。

下一篇是model selection與model averaging。