分類

• • 第四章：分 類
• 4.1 分類問題叙述
• 4.2 為什麼線性回歸不可用
• 4.3 Logistic Regression（邏輯回歸）
• 4.4 線性判别分析（LDA）
• 4.5 二次判别分析（QDA）
• 4.6 分類方法的比較
• • 總而言之，選擇壓倒一切，沒有萬能的模型，隻有适合的模型。
• 4.7 具體代碼實作（R語言版本）
• • 4.7.1 R語言版本
• 4.8 結語

第四章：分 類

4.1 分類問題叙述

• 線上性回歸模型中，我們假設相應變量Y是定量的，如：Y是從2001年到2019年中國總人口數量等。但是在很多情況下，響應變量Y 卻是 定性的(qualitative)，例如：人的眼睛顔色是定性變量，取值為藍色，棕色，綠色。定性變量也稱為 分類(categorical) 變量，兩者的統計意義是一樣的。

• 現實生活中分類問題是比較常見的，下面舉幾個例子：

  (1). 某人來到急診室，醫生通過病人的一系列症狀将病人歸類到三類可能的病症中的一類。

  (2). 網上銀行服務需要基于使用者的IP位址，交易記錄等資訊辨識一個線上網上交易是否存在詐騙。

  (3). 一位生物學家對一定數量的患先天疾病的患者與沒有患有先天疾病的人進行DNA排序資訊分析，發現哪些的DNA突變是緻病的，哪些是不緻病的。

• 下面我們将用一個模拟的資料集Default(違約)資料集來闡述分類模型的概念。

  在該問題中，感興趣的是一個人的年收入和月信用卡餘額來預測其違約的狀态資料集如圖4-1所示：（橙色為違約狀态，藍色為未構成違約狀态）

  

  

4.2 為什麼線性回歸不可用

在之前，我們已經用線性回歸解決了很多問題，那我們自然就能想到線性回歸能不能用于解決分類問題呢？那我們來探讨這個極具價值的問題吧：

• 我們先用線性回歸來進行模組化： d e f a u l t = β 0 + β 1 B a l a n c e + β 2 I n c o m e {default = \beta_0 + \beta_1 Balance + \beta_2 Income} default=β0​+β1​Balance+β2​Income，隻要輸入Balance 和 Income 以及default的資料就能用最小二乘法估計出 β 0 , β 1 {\beta_0,\beta_1} β0​,β1​,設定預測的default>0.5就是違約反之不違約，感覺很完美的樣子，但事實真的是這樣嗎？

• 我們仔細來看這個線性回歸模型：

  (1).我們假設有一個窮人Lisa,他的Balance和Income都很小，那麼有可能會導緻default的值為負數，那麼這個負數代表什麼意義呢？

  

  （圖4-2表示線性回歸模型，在這個模型上有些default竟然是負的！）

  (2).當我們的分類變量是多類的時候，以0.5為界限劃分分類就不可用了，那麼我們應該怎麼找到一個界限衡量多分類呢？

  基于以上問題，現在大家是否還覺得線性回歸模型作為一個分類模型是否足夠優秀呢？其實，為了解決以上的問題（1）我們來想想能不能将線性回歸的結果default轉化為區間[0:1]上，讓default轉變成一個違約的機率呢？下面我們來解決這個問題吧。

4.3 Logistic Regression（邏輯回歸）

• 在推導邏輯回歸之前，我們先來認識下一組函數，這組函數具有神奇的作用，可以将是實數軸上的數轉換為[0:1]區間上的機率。

  首先，我們假設我們的線性回歸模型為 Y = β 0 + β 1 X {Y=\beta_0+\beta_1 X} Y=β0​+β1​X，那麼這個函數是如何将線性回歸的結果轉化為機率呢？這個函數就是logistic 函數，具體的形式為 p ( X ) = e β 0 + β 1 X 1 + e β 0 + β 1 X {p(X) = \dfrac{e^{\beta_0 + \beta_1X}}{1+e^{\beta_0 + \beta_1X}}} p(X)=1+eβ0​+β1​Xeβ0​+β1​X​，他的函數圖像如圖4-2：（左邊是線性回歸，右邊是邏輯函數）

• 那麼問題又來了，我們怎麼去估計 p ( X ) = e β 0 + β 1 X 1 + e β 0 + β 1 X {p(X) = \dfrac{e^{\beta_0 + \beta_1X}}{1+e^{\beta_0 + \beta_1X}}} p(X)=1+eβ0​+β1​Xeβ0​+β1​X​ 中的 β 0 , β 1 {\beta_0 , \beta_1} β0​,β1​呢？顯然，極大似然估計是個好辦法，我們先寫出 p ( X ) = e β 0 + β 1 X 1 + e β 0 + β 1 X {p(X) = \dfrac{e^{\beta_0 + \beta_1X}}{1+e^{\beta_0 + \beta_1X}}} p(X)=1+eβ0​+β1​Xeβ0​+β1​X​的似然函數，如果不懂似然函數的可以用檢視機率統計，這裡直接用結論：我們寫出他的似然函數 l ( β 0 , β 1 ) = ∏ i : y i = 1 p ( x i ) ∏ i : y i ‘ = 0 ( 1 − p ( x i ) ) {\mathcal{l(\beta_0,\beta_1) = {\prod\limits_{i:y_{i} =1}p(x_i)}} {\prod\limits_{i:y_{i ^ `}=0}(1-p(x_i))}} l(β0​,β1​)=i:yi​=1∏​p(xi​)i:yi‘​=0∏​(1−p(xi​))，對這個似然函數求對數，對 β 0 , β 1 {\beta_0,\beta_1} β0​,β1​求偏導，求極值就可以得出系數啦，不過現實中還是使用相關的庫求吧，手算太麻煩了，結果如下圖所示：（這裡沒有将income也寫進去）

  

• 到現在為止，我們已經推導出這個邏輯回歸模型了，并且已經把所有的系數都估計出來啦，離問題解決隻剩下最後一步了，那就是預測。其實，預測非常easy，隻要将資料扔進去我們剛剛求出來的模型就好啦，下面我們來看看吧：

  加入Lisa的Balance=1000美元時，我們剛剛建立的模型是： p ^ ( X ) = e β 0 ^ + β 1 ^ X 1 + e β 0 ^ + β 1 ^ X {\hat{p}(X) = \dfrac{e^{\hat{\beta_0} + \hat{\beta_1 }X}}{1+e^{\hat{\beta_0} + \hat{\beta_1 }X}}} p^​(X)=1+eβ0​^​+β1​^​Xeβ0​^​+β1​^​X​，我們把balance=1000代入X中，得出 p ^ ( X ) = 0.00576 {\hat{p}(X) = 0.00576} p^​(X)=0.00576，違約的機率小于1%，可見她不會違約咯。

• 基于以上問題我們已經成功地将線性回歸模型過渡到分類模型，讓模型更加适合分類，我們成功地解決了問題(1)，但是還有點小瑕疵，那就是多個自變量的時候我們沒有說明白，其實非常簡單，就是将一進制線性回歸變成多元線性回歸輸入我們的邏輯函數， p ( X ) = e β 0 + β 1 X 1 + β 2 X 2 + β 3 X 3 + . . . . 1 + e β 0 + β 1 X 1 + β 2 X 2 + β 3 X 3 + . . . . {p(X) = \dfrac{e^{\beta_0 + \beta_1X_1+\beta_2X_2+\beta_3X_3+....}}{1+e^{\beta_0 + \beta_1X_1+\beta_2X_2+\beta_3X_3+....}}} p(X)=1+eβ0​+β1​X1​+β2​X2​+β3​X3​+....eβ0​+β1​X1​+β2​X2​+β3​X3​+....​，就是這麼簡單咯。
• 值得注意的一點是：建立完模型記得看P-value的大小，P-value < 0.05才可以使用這個模型，如果不滿足隻能選擇其他變量建立模型了。

  

• 對于問題(2),我們值得注意的是，邏輯回歸在實際中不太用于多分類問題，因為實際效果不是很好，是以我們必須借助其他模型來解決這個問題，那讓我們來解決這個遺留下來的問題吧。

4.4 線性判别分析（LDA）

• 在讨論如何解決多分類問題之前，我們先來說說貝葉斯的那些事吧。在機率統計的領域裡有一條神奇的公式叫貝葉斯定理，具體的形式是： P ( Y = k ∣ X = x ) = π k f k ( x ) ∑ l = 1 K π l f l ( x ) {P(Y=k|X=x) = \dfrac{{\pi}_kf_k(x)}{\sum\limits_{l=1}^K{\pi}_lf_l(x)}} P(Y=k∣X=x)=l=1∑K​πl​fl​(x)πk​fk​(x)​ ，我們先不要被公式的符号吓到，我們先來看看符号具體代表什麼意思。我們假設觀測有 K {K} K類， π k {\pi_k} πk​為随機選擇的觀測來自第 k {k} k類的 先驗機率，也就是樣本裡面第 k {k} k類的樣本個數除以總樣本的個數： π k = n k n {\pi_k = \dfrac{n_k}{n}} πk​=nnk​​。再來 f k ( x ) = P ( X = x ∣ Y = k ) {f_k(x) =P(X=x|Y=k)} fk​(x)=P(X=x∣Y=k)，表示第 k {k} k類觀測的X的密度函數，說的直白一點就是在 Y = k {Y=k} Y=k的樣本裡 X = x {X=x} X=x的樣本個數，即 f k ( x ) = P ( X = x ∣ Y = k ) = n ( X = x , Y = k ) n ( Y = k ) {f_k(x) = P(X=x|Y=k) = \dfrac{n_{(X=x,Y=k)}}{n_{(Y=k)}}} fk​(x)=P(X=x∣Y=k)=n(Y=k)​n(X=x,Y=k)​​，最後， ∑ l = 1 K π l f l ( x ) = P ( X = x ) = n ( X = x ) n {\sum\limits_{l=1}^K{\pi}_lf_l(x)}=P(X=x)=\dfrac{n_{(X=x)}}{n} l=1∑K​πl​fl​(x)=P(X=x)=nn(X=x)​​，也就是樣本中 X = x {X=x} X=x的機率。
• 在讨論貝葉斯定理後，我們回到分類問題，這個定理跟我們的分類問題有什麼關聯呢？沒錯，這個公式 P ( Y = k ∣ X = x ) = π k f k ( x ) ∑ l = 1 K π l f l ( x ) {P(Y=k|X=x) = \dfrac{{\pi}_kf_k(x)}{\sum\limits_{l=1}^K{\pi}_lf_l(x)}} P(Y=k∣X=x)=l=1∑K​πl​fl​(x)πk​fk​(x)​給出了給定樣本條件下， Y = k {Y=k} Y=k這個類别下的機率，這給分類問題提供了一條思路，那就是計算這個 P ( Y = k ∣ X = x ) {P(Y=k|X=x)} P(Y=k∣X=x)，而且我們的邏輯回歸就是這麼幹的，但是在 P ( Y = k ∣ X = x ) = π k f k ( x ) ∑ l = 1 K π l f l ( x ) {P(Y=k|X=x) = \dfrac{{\pi}_kf_k(x)}{\sum\limits_{l=1}^K{\pi}_lf_l(x)}} P(Y=k∣X=x)=l=1∑K​πl​fl​(x)πk​fk​(x)​這個公式中，分母 ∑ l = 1 K π l f l ( x ) = P ( X = x ) {{\sum\limits_{l=1}^K{\pi}_lf_l(x)} = P(X=x)} l=1∑K​πl​fl​(x)=P(X=x)當樣本給定的時候是一個與分類 k {k} k無關的常數,是以我們的問題可以簡化為隻需要計算分子 π k f k ( x ) {{\pi}_kf_k(x)} πk​fk​(x),進而比較哪個類别的機率最大就知道屬于哪個類别了，是以我們的分類思路就出來啦，這個思路不同于邏輯回歸，邏輯回歸需要計算具體的 P ( Y = k ∣ X = x ) {P(Y=k|X=x)} P(Y=k∣X=x)機率值，而我們現在的思路是通過貝葉斯定理計算貝葉斯定理的分子，比較分子最大的那個類别為最終類别。

• 在我們推導複雜算法之前，我們先推導下簡單的當自變量個數隻有一個的模型，即 p = 1 {p=1} p=1的簡單模型。我們記 P ( Y = k ∣ X = x ) = π k f k ( x ) ∑ l = 1 K π l f l ( x ) {P(Y=k|X=x) = \dfrac{{\pi}_kf_k(x)}{\sum\limits_{l=1}^K{\pi}_lf_l(x)}} P(Y=k∣X=x)=l=1∑K​πl​fl​(x)πk​fk​(x)​ 的分子為 g k ( x ) = π k f k ( x ) {g_k(x) = {\pi}_kf_k(x)} gk​(x)=πk​fk​(x)。在這裡，我們做個模型假設：假設 f k ( x ) {f_k(x) } fk​(x)服從正态分布，即 f k ( x ) ∼ N ( μ , σ k 2 ) {f_k(x) \sim N(\mu,\sigma_k^2)} fk​(x)∼N(μ,σk2​)，而且每個 σ k 2 = σ 2 {\sigma_k^2 = \sigma^2} σk2​=σ2，同方差假設。是以 f k ( x ) = 1 2 π σ k e − 1 2 σ 2 ( x − μ k ) 2 {f_k(x) = \dfrac{1}{\sqrt{2\pi}\sigma_k}e^{-\dfrac{1}{2\sigma^2}(x-\mu_k)^2}} fk​(x)=2π

  ​σk​1​e−2σ21​(x−μk​)2，最終我們的 g k ( x ) = π k 1 2 π σ k e − 1 2 σ 2 ( x − μ k ) 2 {g_k(x) = \pi_k\dfrac{1}{\sqrt{2\pi}\sigma_k}e^{-\dfrac{1}{2\sigma^2}(x-\mu_k)^2}} gk​(x)=πk​2π

  ​σk​1​e−2σ21​(x−μk​)2,終于算出來啦。這個式子不是很好計算，我們對 g k ( x ) {g_k(x)} gk​(x)取個對數，令 δ k ( x ) = l n ( g k ( x ) ) = l n π k + μ σ 2 x − μ 2 2 σ 2 {\delta_k(x) = ln(g_k(x))=ln\pi_k+\dfrac{\mu}{\sigma^2}x-\dfrac{\mu^2}{2\sigma^2}} δk​(x)=ln(gk​(x))=lnπk​+σ2μ​x−2σ2μ2​，到這裡我們的模型建立模型，我們隻需要把位置的 μ k {\mu_k} μk​與 σ 2 {\sigma^2} σ2估計出來就好了。 μ ^ k = 1 n k ∑ i : y i = k x i {\hat{\mu}_k =\dfrac{1}{n_k}\sum\limits_{i:y_i=k}x_i} μ^​k​=nk​1​i:yi​=k∑​xi​，也就是當 y = k {y=k} y=k這一類中 x {x} x的平均值； σ ^ 2 = 1 n − K ∑ k = 1 K ∑ i : y i = k ( x i − μ ^ k ) 2 {\hat{\sigma}^2 =\dfrac{1}{n-K}\sum\limits_{k=1}^K\sum\limits_{i:y_i=k}(x_i-\hat{\mu}_k)^2 } σ^2=n−K1​k=1∑K​i:yi​=k∑​(xi​−μ^​k​)2，說白了就是計算每一類的方差，再求平均值。總結下上面的公式就是：

  { δ k ( x ) = l n ( g k ( x ) ) = l n π k + μ σ 2 x − μ 2 2 σ 2 μ ^ k = 1 n k ∑ i : y i = k x i σ ^ 2 = 1 n − K ∑ k = 1 K ∑ i : y i = k ( x i − μ ^ k ) 2 {\begin{cases}\delta_k(x) = ln(g_k(x))=ln\pi_k+\dfrac{\mu}{\sigma^2}x-\dfrac{\mu^2}{2\sigma^2}\\{\hat{\mu}_k =\dfrac{1}{n_k}\sum\limits_{i:y_i=k}x_i}\\{\hat{\sigma}^2 =\dfrac{1}{n-K}\sum\limits_{k=1}^K\sum\limits_{i:y_i=k}(x_i-\hat{\mu}_k)^2}\end{cases}} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​δk​(x)=ln(gk​(x))=lnπk​+σ2μ​x−2σ2μ2​μ^​k​=nk​1​i:yi​=k∑​xi​σ^2=n−K1​k=1∑K​i:yi​=k∑​(xi​−μ^​k​)2​

  至此，我們的模型就建立完成了，我們隻需要代入資料求出 δ k ( x ) {\delta_k(x)} δk​(x)，哪個 k {k} k對應的 δ k ( x ) {\delta_k(x)} δk​(x)大，就是哪一類。

  （下圖虛線是線性判别分析的決策邊界，正态曲線哪邊高樣本就是哪一類）

  

• 就像邏輯回歸一樣，我們推到出了一個自變量的簡單模型，就要泛化為多個自變量的線性判别分析了，即 p > 1 {p>1} p>1。其實原理一樣的，隻是将一進制正态分布擴充為多元正态分布：

  f k ( x ) = 1 2 ( π ) p 2 ∣ Σ ∣ 1 2 e [ − 1 2 ( x − μ k ) T Σ ( x − μ k ) ] {f_k(x)=\dfrac{1}{2(\pi)^{\tfrac{p}{2}}|\Sigma|^\tfrac{1}{2}}e^{[-\tfrac{1}{2}(x-\mu_k)^T\Sigma(x-\mu_k)]}} fk​(x)=2(π)2p​∣Σ∣21​1​e[−21​(x−μk​)TΣ(x−μk​)]

  μ k ^ = ( μ k 1 , μ k 2 , . . . . . . , μ k p ) , Σ ^ = 1 p − 1 ∑ j = 1 p ( x j − x ‾ ) ( x j − x ‾ ) T {\hat{\mu_k}=(\mu_{k1},\mu_{k2},......,\mu_{kp}) , \hat{\Sigma}=\dfrac{1}{p-1}\sum\limits_{j=1}^p(x_j-\overline{x})(x_j-\overline{x})^T} μk​^​=(μk1​,μk2​,......,μkp​),Σ^=p−11​j=1∑p​(xj​−x)(xj​−x)T

  δ k ( x ) = l n ( π k f k ( x ) ) = l n ( π k ) − ( p 2 l n ( 2 π ) + 1 2 l n ( ∣ Σ ∣ ) ) − 1 2 ( x − μ k ) T Σ − 1 ( x − μ k ) = x T Σ ^ μ ^ k − 1 2 μ ^ k T Σ ^ − 1 μ ^ k + l n π ^ k {\delta_k(x) = ln(\pi_kf_k(x))=ln(\pi_k)-(\dfrac{p}{2}ln(2\pi)+\dfrac{1}{2}ln(|\Sigma|))-\dfrac{1}{2}(x-\mu_k)^T\Sigma^-1(x-\mu_k)=x^T\hat{\Sigma}\hat{\mu}_k-\dfrac{1}{2}\hat{\mu}_k^T\hat{\Sigma}^{-1}\hat{\mu}_k+ln\hat{\pi}_k} δk​(x)=ln(πk​fk​(x))=ln(πk​)−(2p​ln(2π)+21​ln(∣Σ∣))−21​(x−μk​)TΣ−1(x−μk​)=xTΣ^μ^​k​−21​μ^​kT​Σ^−1μ^​k​+lnπ^k​

• 從上面的分析，我們至少能完整解決多分類問題了，那我們又開始想一個問題了，能不能優化這個算法呢？線性判别分析究竟存在什麼問題呢？

4.5 二次判别分析（QDA）

• 為了優化線性判别分析，我們從線性判别分析的模型假設入手，我們假設 f k ( x ) {f_k(x)} fk​(x)是同方差的正态分布，然而現實中真的完全都是同方差嗎？例如同是藍色眼睛的人身高和體重和黑眼睛的身高體重的方差一定一樣嗎？為了打破這個假設，我們必須假設 f k ( x ) {f_k(x)} fk​(x)的方差各自不同，即我們假設 Σ 1 ≠ Σ 2 ≠ . . . . . . ≠ Σ K {\Sigma_1\neq \Sigma_2\neq......\neq\Sigma_K} Σ1​​=Σ2​​=......​=ΣK​,那判别函數就變成 δ k ( x ) = − 1 2 x T Σ k − 1 x + x T Σ k − 1 μ k − 1 2 μ k T Σ k − 1 μ k + l n π k {\delta_k(x) = -\dfrac{1}{2}x^T\Sigma^{-1}_kx+x^T\Sigma^{-1}_k\mu_k-\dfrac{1}{2}\mu^{T}_k\Sigma^{-1}_k\mu_k+ln\pi_k} δk​(x)=−21​xTΣk−1​x+xTΣk−1​μk​−21​μkT​Σk−1​μk​+lnπk​，我們仔細觀察下這個表達式，就會知道 當不同 的 C o v ( x i , x j ) = 0 {Cov(x_i,x_j)=0} Cov(xi​,xj​)=0時，剛好QDA模型就是樸素貝葉斯模型，是以樸素貝葉斯是簡化版的QDA。

  （下圖是QDA的綠色決策邊界，可見QDA模型是非線性的二次函數決策邊界）

  

  

4.6 分類方法的比較

• 在前面的探讨中，我們已經完美地解決了某些分類問題，并且我們從兩個大的角度：計算 P ( Y = k ∣ X = x ) {P(Y=k|X=x)} P(Y=k∣X=x)還是比較 δ k ( x ) {\delta_k(x)} δk​(x)的大小來決定某個樣本的類别，也建立了三個分類模型：邏輯回歸，線性判别分析LDA和二次判别分析QDA來解決分類問題，同時之前我們也探讨過了K近鄰法（KNN）也可以解決某些分類問題，那麼在這些方法之間有什麼優劣之處呢？

• 邏輯回歸 VS 線性判别分析LDA

  在回答這個問題之前，我們來重新對這兩個模型的公式做些比較。首先是LDA，我們假設類别 K = 2 {K=2} K=2，并且假設每種類别分别為 K = 1 ， K = 2 {K=1，K=2} K=1，K=2，設 P 1 ( X ) , P 2 ( X ) = 1 − P 1 ( X ) {P_1(X),P_2(X)=1-P_1(X)} P1​(X),P2​(X)=1−P1​(X)分别為 K = 1 ， K = 2 {K=1，K=2} K=1，K=2的機率。對LDA模型下取對數發生比 l o g ( P 1 ( X ) 1 − P 1 ( X ) ) = l o g ( P 1 ( X ) P 2 ( X ) ) = c 0 + c 1 X {log(\dfrac{P_1(X)}{1-P_1(X)})=log(\dfrac{P_1(X)}{P_2(X)})=c_0+c_1X} log(1−P1​(X)P1​(X)​)=log(P2​(X)P1​(X)​)=c0​+c1​X，其中 c 0 , c 1 {c_0,c_1} c0​,c1​是 μ 0 , μ 1 , σ 2 {\mu_0,\mu_1,\sigma^2} μ0​,μ1​,σ2的函數。

  再來看邏輯回歸， l o g ( P 1 P 2 ) = β 0 + β 1 X {log(\dfrac{P_1}{P_2})=\beta_0+\beta_1X} log(P2​P1​​)=β0​+β1​X，LDA與邏輯回歸的對數發生比都是關于 X {X} X的線性函數，是以他們的決策邊界都是線性的，但由于邏輯回歸使用極大似然估計計算出最優的 β 0 , β 1 {\beta_0,\beta_1} β0​,β1​，而LDA是用正态分布的均值和方差計算出來的，他們的拟合過程有差異，但是他們的結果應該是相近的，但并非是必然的。LDA假設觀測服從每一類協方差矩陣都相同的正态分布（均值不一定相同），是以當樣本對這個假設近似滿足的時候，LDA應該會比邏輯回歸的效果會更好，如果假設不滿足，則邏輯回歸的效果比LDA更好。

• KNN VS 邏輯回歸與LDA

  當我們談論LDA時，我們應該清楚KNN的分類原理跟我們剛剛談論的完全不同，對觀測 X = x {X=x} X=x進行預測時，先找出該點最近的 K {K} K個點，找出K個點中屬于類别數最多的點作為該點的分類，說白了就是“少數服從多數”，是以KNN是一個徹底的非參數方法，也就是他沒有假設決策邊界的形狀。是以，當決策邊界高度非線性時，KNN會比邏輯回歸和LDA效果更好，反之。

• 二次判别分析QDA VS KNN,邏輯回歸，LDA

  QDA是非參數方法KNN與線性方法LDA，邏輯回歸的一種折中方法，因為QDA假設得到的是一個二次決策邊界，他比線性決策邊界的模型适用性更廣，雖然不如KNN方法的光滑程度（模型複雜度）高，但是QDA在固定訓練資料量的問題上一般比KNN有更好的效果，因為QDA對決策邊界的形狀做了一定的假設。

  （下圖可看到在不通類型的資料下，各類模型的錯誤率不一樣。）

  

• 最後做個總結，當真實決策邊界是線性的LDA和邏輯回歸會比較好，當真實決策邊界是一般非線性的QDA會比較好，當真實決策邊界是更加複雜的，KNN會比較合适。

總而言之，選擇壓倒一切，沒有萬能的模型，隻有适合的模型。

4.7 具體代碼實作（R語言版本）

4.7.1 R語言版本

為什麼使用R而不用python呢？因為這幾個我們的模型都是傳統的統計模型，統計模型在R這種專業的統計語言裡比較全面，是以在面對統計模型時R是個十分不錯的選擇。

• 1.資料的初步探索：

  我們将會用到兩個資料: Smarket資料和Purchase資料。這兩個資料都在R包ISLR中，是以需要我們先用函數require()或者library()加載ISLR(如果電腦中沒有安裝ISLR，請先安裝ISLR，然後再加載。首先，請先使用函數class(), names(), dim(), head(), summary()等初步探索該資料，可以對資料有一個直覺的認識。例如：

require(ISLR)   # 加載ISIR包
names(Smarket)  # 檢視Smarket的變量名
dim(Smarket)    # 檢視Smarket資料集的次元
cor(Smarket[, -9])   # 檢視各變量的相關性
           

如果大家對股票投資有一些了解，可以看出上面的結果和一般的預期是一緻的，大部分變量之間的相關性都是非常小的。如果發現某些變量和當天的收益率相關性很強，投資者就可以利用這種相關性投資股票，獲利豐厚了。現實是，在股票上賺錢是一件很困難的事情，是以也不太可能這麼容易就可以找到和未來收益相關的變量。仔細看上面的結果，也可以發現“Year”和“Volume”之間的相關關系是“0.539”, 一個比較大的相關關系。這是為什麼呢？

require("RColorBrewer")    # 設定顔色
miscolores = brewer.pal(8,"Set2")[1:8]  # 設定顔色
par(mar=c(5,5,0.1,0.1))   # 設定顔色
plot(Smarket[, "Volume"], col = miscolores[1], xlab = "Index", ylab = "Volume")  # 作出Volume散點圖
           

可以發現，“Volume”總的趨勢是在上升，說明該股票的交易越來越活躍。是以，“Year”和“Volume”可以計算出一個比較大的相關關系。

直覺來看，Sarmket資料中有兩個變量可以作為因變量，“Today”和“Direction”。“Today”表示股票當天的收益率，是定量變量。“Direction”表示股票當天是漲還是跌，是定性變量。是以，當“Today”是非負的，“Direction”等于“Up”；反之，當“Today”是負的，“Direction”等于“Down”。

attach(Smarket) # 引用Smarket資料，以便在後面不需要再同一個代碼塊内重複輸入資料名字
# 把Today>0的變量direction設為up，反之為down
my_Direction <- ifelse(Today >= 0, "Up", "Down")
sum(my_Direction == Direction)
as.factor(my_Direction[1:5])
Direction[1:5]
detach() # 解除資料約定
           

當“Today”作為因變量時，我們需要建立的是回歸模型；當“Direction”作為因變量時，我們需要建立分類模型。在這次的實驗中，我們要學習的是分類方法的模組化，是以把“Direction”作為因變量。

• 2.Logistic回歸模型

  2.1 建立邏輯回歸模型

  在R中，可以使用函數glm()建立logistic回歸模型，函數glm()和函數lm()的使用方法有很多一樣的地方。特别需要注意的是，函數glm()需要使用參數family = binomial來表示我們要建立的模型是logristic回歸模型。因為, “glm”表示generalized linear model(廣義線性模型), 函數glm()可以拟合這個系列的模型。同學們可以在R中輸入?glm檢視其它的參數family的選擇。

glm.fit <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Smarket, family = binomial) # 建立邏輯回歸模型
glm.fit
           

我們建立的logistic模型，以Direction為因變量，以Lag1, Lag2, Lag3, Lag4, Lag5, Volume為自變量。輸入glm.fit可以看到logistic模型的拟合系數。和函數lm()的結果一樣， 我們也可以用names()可以看到glm.fit包含哪些資訊：

我們也可以用函數summary()看到更多的資訊：

函數summary()的結果主要包含了系數，系數的标準差，統計量，p值。從上面的結果，可以看到所有的自變量的p-value都比較大，說明這些自變量對Direction沒有顯著的作用。在這些自變量中，Lag1的p-value最小，雖然它的p-value還是比較大，等于0.15。Lag1的系數的值是負的，說明了如果前一天股票是下跌的，那麼該股票今天上漲的可能性比較大；如果前一天股票是上漲的，那麼該股票今天下跌的可能性比較大。

這模型glm.fit中，我們知道R在計算的過程中把Direction的Up或者Down變成了“1”。從上述的模型中我們無法得到這個資訊。這時我們可以使用函數contrasts()檢視R是怎麼把定性變量變成傀儡變量的：

在Smarkdet資料中，R把Up記為1， Down記為0。 函數predict()可以得到得到機率P(Y=1|X)。在這個例子中，函數predict()得到的是股票上漲的機率。

glm.probs <- predict(glm.fit, type = "response")
length(glm.probs)  # 檢視結果的長度
           

如果大家在R console中輸入?predict.glm檢視函數predict()的幫助文檔，可以看到第一個參數是object, 該object是函數glm()的傳回值，在我們這裡也就是“glm.fit”； 第二個參數是newdata, 指的是我們要預測的自變量的值組成的資料框，該參數的預設值是訓練資料，即函數predict()的傳回值是拟合值；第三個參數是type, 有三個選擇，“link”, “response”, “terms”, 一般我們選擇“response”, 表示我們要得到的結果是機率P(Y=1|X)。在這裡，我們沒有改變第二個參數newdata的值，是以，函數predict()預設使用訓練資料作為參數newdata的輸入，是以我們可以看到訓練資料中1250個觀測值的拟合機率P(Y=1|X)。結果“glm.probs”的第一個數字表示第一個觀測值的預測機率為0.5070841。

如果我們已經知道昨天的成交量“Volume = 1.5”, 并且“Lag1 = 1, Lag2 = 1, Lag3 = 1, Lag4 = 1, Lag5 = 1”, 我們想預測今天的股票是漲還是跌，可以使用如下的方式：

glm.today <- predict(glm.fit, newdata = data.frame(Lag1 = 1, Lag2 = 1, Lag3 = 1, Lag4 = 1, Lag5 = 1, Volume = 1.5), type = "response")
glm.today
           

通過使用函數predict()，我們已經得到了機率值P(Y=1|X)。如果我們使用如下規則判斷股票漲跌：如果機率值P(Y=1|X)大于等于0.5，那麼我們判斷今天的股票是漲,“Up”;如果如果機率值P(Y=1|X)小于0.5，那麼我們判斷今天的股票是跌,“Down”。可以使用下面兩行代碼實作上面的規則:

glm.pred <- rep("Down", 1250)
glm.pred[glm.probs >= 0.5] <- "Up"
glm.pred[1:6]
           

第一行代碼表示産生一個向量，該向量的元素全部是“Down”, 向量的長度是1250。第二行代碼glm.probs >= 0.5判斷哪些觀測點的“Up”的機率大于等于0.5。然後把glm.pred中對應的元素指派為“Up”。然後我們可以使用函數table()得到一個confusion matrix:

attach(Smarket)
table(glm.pred, Direction)
detach()
(507 + 145) / 1250
mean(glm.pred == Smarket$Direction)
           

函數table()輸出的結果就是confusion matrix了。在上面的代碼中，函數table()的第一個參數是預測的股票漲跌，第二個參數是實際的股票漲跌。輸出的結果confusion matrix，清單示預測的結果，行表示實際的結果。是以，confusion matrix中，145表示實際股票是“Down”,預測也為“Down”的觀測點個數是145; 141表示實際股票是“Up”, 但是預測為“Down”的觀測點個數是141；457表示實際股票是“Down”, 但是預測為“Up”的觀測點個數是457；507表示實際股票是“Up”,預測也為“Up”的觀測點個數507。是以正确的預測比例是0.5216。

在股票的預測中，僅僅從幾天内的收益率與成交量的資訊，用logistic 回歸模型得到52%的成績應該已經算不錯的成績了。但是，值得注意的是現在得到的這個正确率是訓練資料的拟合正确率。該模型的訓練資料的錯誤率是100%−52%=48%。我們說過，訓練誤差(training error)一般情況下是過于樂觀的，也就是說，訓練誤差通常是小于真正的預測誤差的。是以，現在得到的這個訓練正确率有可能是大于真實的正确率的。

2.2Logistic回歸模型的預測誤差

為了得到更加準确的預測正确率的估計，一個自然的想法是，隻用一部分資料建立Logistic回歸模型，然後用剩下的資料計算模型的預測誤差。在這個資料中，我們可以使用2001-2004的資料建立一個Logistic回歸模型，然後再用建立的模型預測2005的股票漲跌，把預測的股票漲跌與真實的股票漲跌相比較，我們就可以得到模型的預測正确率或者錯誤率。

train <- (Smarket$Year < 2005)
train[1:10]
(!train)[1:10]
summary(train)
           

Smarket_2001TO2004 <- Smarket[train, ]
head(Smarket_2001TO2004)
tail(Smarket_2001TO2004)
           

Smarket_2005 <- Smarket[!train, ]
head(Smarket_2005)
           

train是一個布爾向量(Boolean vector), train的長度是1250，隻包含FALSE或者“TRUE”。在這裡，如果資料的Year小于2005，那麼對應位置的train的值指派為“TRUE”, 如果資料的Year大于2005，那麼對應位置的train的值指派為“FALSE”。在R中， !表示相反的，即!TRUE表示FALSE, !FALSE表示TRUE。是以, !train表示與train長度相同的向量，隻是train向量中的TRUE，在!train中是FALSE; train向量中的FALSE，在!train中是TRUE。 表達式 Smarket[train, ]表示取資料框 Smarket中train為TRUE所對應的行。是以，Smarket_2001TO2004表示2001-2004的股票資料， Smarket_2005表示2005的股票資料。

這時我們可以使用2005前的股票資料建立Logistic回歸模型, 然後再預測2005的股票漲跌。隻需要在函數glm()中加入參數subset:

glm.fit <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Smarket, family = binomial, subset = train)
summary(glm.fit)
           

然後預測2005的股票漲跌:

glm.probs <- predict(glm.fit, newdata = Smarket_2005, type = "response")
length(glm.probs)
           

glm.preds <- rep("Down", length(glm.probs))
glm.preds[glm.probs >= 0.5] <- "Up"
glm.probs[1:6]
glm.preds[1:6]
           

現在，我們可以用函數table()得到confusion matrix:

table(glm.preds, Smarket_2005$Direction)
mean(glm.preds == Smarket_2005$Direction)
mean(glm.preds != Smarket_2005$Direction)
           

在R中，符号== 表示判斷 ==兩邊的向量或者矩陣對應的元素是否相等；符号!=表示判斷!=兩邊的向量或者矩陣對應的元素是否不相等。是以, mean(glm.preds == Smarket_2005 D i r e c t i o n ) 表 示 預 測 正 确 率 ； m e a n ( g l m . p r e d s ! = S m a r k e t 2 005 Direction)表示預測正确率；mean(glm.preds != Smarket_2005 Direction)表示預測正确率；mean(glm.preds!=Smarket2​005Direction)表示預測錯誤率。

現在，預測正确率隻有48%, 低于扔硬币的準确率！說明我們建立的Logistic回歸模型不僅不能幫助我們預測股票走勢，正确率反而低于随機的猜測。

有什麼簡單的辦法可以幫助我們提高預測的準确率呢？我們認真看看上面summary(glm.fit)的結果，可以看到，所有的自變量的系數都不是顯著的不等于0的，即，所有自變量的p值都是大于0.05的。隻是Lag1和Lag2的p值相對其他的小一些。如果我們隻選取Lag1和Lag2作為自變量建立Logistic回歸模型效果會不會更好一些呢？

glm.sub <- glm(Direction ~ Lag1 + Lag2, data = Smarket_2001TO2004, family = binomial)
glm.sub.probs <- predict(glm.sub, newdata = Smarket_2005, type = "response")
glm.sub.preds <- rep("Down", length(glm.sub.probs))
glm.sub.preds[glm.sub.probs >= 0.5] <- "Up"
table(glm.sub.preds, Smarket_2005$Direction)
mean(glm.sub.preds == Smarket_2005$Direction)
           

這時我們的模型表現的比之前的好很多了，正确率可以達到56%。從confusion matrix, 我們也可以看到，當模型預測股票跌時，模型的正确率隻有35/(35+35)=50%; 當模型預測股票漲時，模型的正确率有106/(106+76)=58%。是以，我們可以采用如下的政策投資股票：模型預測為跌時，我們不做任何的投資；模型預測為漲時，我們可以買進股票。這個政策的正确率可以達到58%!

• 3. LDA與QDA

  3.1 LDA線性判别分析

  在這一小節中，我們将學習使用函數lda()在R中實作LDA。函數lda()在R的包MASS中，并且函數lda()的使用方法和函數lm()或者glm()是非常相似的。

require(MASS)
lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket_2001TO2004)
lda.fit
           

lda.fit包含了很多有用的資訊，比如，每一類的先驗機率。上面結果中，0.491984表示在資料Smarket_2001TO2004中，有49.1984%天該股票是下跌的。

與函數lm()的結果類似，我們也可以使用函數predict()得到預測值:

lda.pred <- predict(lda.fit, Smarket_2005)
names(lda.pred)
           

結果lda.pred包含了最終的預測值，class和posterior。posterior表示預測為“Down”或者“Up”的機率。class是使用門檻值0.5對posterior分類的結果，即，當posterior中，表示“Up”的機率大于0.5時，class = “Up”， 反之，class = “Down”。

lda.pred$class[1:10]
head(lda.pred$posterior)
post.pred <- rep("Down", nrow(Smarket_2005))
post.pred[lda.pred$posterior[, "Up"] > 0.5] <- "Up"
sum(post.pred == lda.pred$class)
nrow(Smarket_2005)
           

這裡，同樣的我們可以用函數table()計算confusion matrix:

table(lda.pred$class, Smarket_2005$Direction)
mean(lda.pred$class == Smarket_2005$Direction)
           

從上面結果可以看到，在資料Smarket中, 使用logistic回歸模型和lda的效果是非常接近的。(印證了我們剛剛的理論推導)

3.2 QDA二次判别分析

在R中，QDA可以使用函數qda()實作，函數qda()同樣包含在R包MASS中。而且使用的方法和函數lda()基本一樣。

qda.fit <- qda(Direction ~ Lag1 + Lag2, data = Smarket_2001TO2004)
qda.fit
           

和函數lda()的結果相似，函數qda()的結果也包含了不同類的先驗機率和每一個類下，自變量的平均值。可以使用函數predict()得到預測值:

qda.pred <- predict(qda.fit, Smarket_2005)
names(qda.pred)
table(qda.pred$class, Smarket_2005$Direction)
mean(qda.pred$class == Smarket_2005$Direction)
           

在這個資料中，QDA預測的正确率可以達到60%。這在股票的預測中是一個非常不錯的結果。QDA的表現比LDA的表現好很多，說明了LDA假設不同類中的協方差矩陣相等在這個例子中可能是不成立的。

• 4. KNN

  在R中，可以使用包class中的函數knn()實作。函數knn()的文法和我們學過的lm(), lda(), qda()等有所不同。lm(), lda(), qda()等函數做預測都有兩步：第一步用這些函數建立模型，第二步用函數predict()得到預測。而函數knn()需要輸入4個參數:

  (1).一個矩陣或者資料框，包含了訓練資料的自變量的資訊

  (2).一個矩陣或者資料框，包含了測試資料的自變量的資訊

  (3).一個向量，包含了訓練資料的因變量

  (4).一個常數K, 最近的觀測點的個數

首先對資料做一些處理:

train.x <- Smarket_2001TO2004[, c("Lag1", "Lag2")]
test.x <- Smarket_2005[, c("Lag1", "Lag2")]
train.y <- Smarket_2001TO2004[, "Direction"]
           

然後再把這些資料帶入到函數knn()中:

require(class)  # 如果提示沒有包“class”， 要先安裝包“class”哦
knn.pred <- knn(train.x, test.x, train.y, k = 1)
table(knn.pred, Smarket_2005$Direction)
mean(knn.pred == Smarket_2005$Direction)
           

上面的KNN模型，我們使用了參數k=1, 得到的預測正确率是50%。說明模型的效果等同于随機的猜測。我們在課堂上結果，如果k的值很小，模型容易過拟合，最終的預測效果會變差。是以，我們可以嘗試更大的k。

knn.pred <- knn(train.x, test.x, train.y, k = 3)
table(knn.pred, Smarket_2005$Direction)
mean(knn.pred == Smarket_2005$Direction)
knn.pred <- knn(train.x, test.x, train.y, k = 5)
table(knn.pred, Smarket_2005$Direction)
mean(knn.pred == Smarket_2005$Direction)
           

我們可以看到，當k=3時，模型的預測正确率有所提高，但是當k=5時，模型的預測正确率反而降低了。

4.8 結語

到現在為止，我們已經掌握了分類的基本模型，這些模型在我們面對一般的分類問題的時候已經綽綽有餘了，下一步我們将要一起探索幾個問題：

（1）既然我們的模型那麼多，我們應該怎麼在面對資料的時候選出最優秀的模型呢？是邏輯回歸？LDA?QDA?KNN?

（2）當我們標明模型後怎麼選擇合适的變量使我們的模型最優呢？

參考文獻

[1]加雷斯.詹姆斯，丹妮拉.威騰 .統計學習導論[M].北京:機械工業出版社,2015.6.

[2]李航.統計學習方法[M].北京:清華大學出版社,2012.3.

[3]傑克.萬托布拉斯.Python資料科學手冊[M].北京:人民郵電出版社,2018.2.

[4]富朗索瓦.肖萊.Python深度學習[M].北京:人民郵電出版社,2018.8