《信用風險評分卡研究》中最大似然估計分析表的解讀

04-26

眾所周知，Logistic回歸在信用評分卡開發過程中起到了核心的作用，其回歸的結果可以轉換成一個匯總表，形成所謂的標準評分卡，除此之外在《信用風險評分卡研究》一書中還提出了大量對於模型評判的指標，驗證模型的效果，可以說這些指標不僅保證了模型的解釋性也為模型提供了良好的統計學基礎，今天我想談一談的就是模型評判的一些重要指標—最大似然估計分析表，我們先看一下書中第七章出現的兩張表：

了解cost方程&梯度下降

不知道大家有沒有和我有著一樣的疑惑，那就是這些Estimate，Standard Error，Chi-Square，Pr>ChiSq這些參數到底是什麼以及這些參數的怎麼算出來的，他們的含義又都是什麼。今天我就將從一個初學者的角度來解讀一下logistic回歸中這些指標到底怎麼算出來的。

首先我們還是先回顧一下logistic回歸的公式：

其中 $eta_{i}$ 代表每一個x對應的權重，而這個公式最終就是要通過這些x變數和對應的權重去預測最終一行數據輸入後其等於1的概率，通過化簡最終得出：

其中z = $eta^{T}$ x，β和x都是向量，那麼如果當我們最終求出了一個P值以後我們怎麼評判這個結果的好壞呢，這個時候我們一般會使用cost方程，好了第一個問題，什麼是cost方程，讓我們先跳出這個logistic回歸，回到最簡單的線性回歸問題，我們先看一個實際中的簡單例子：

對於如圖所示的一堆散點，最終我們得到了一條紅色的擬合曲線 $h_{ heta}$ = $heta_{0}+ heta_{1}$ x，我們一般會使用如下方式去評判曲線擬合的好壞：

其中J(θ)為擬合的cost方程， $h_{0}(x^{(i)})$ 表示預測值， $y^{(i)}$ 表示實際值，如果擬合的越好那麼J(θ)的值就會越小，我們做機器學習的目的就是要找出合適的 $heta_{0}$ ， $heta_{1}$ 使得這個J(θ)最小，學習過導數的同學都知道，只要讓J(θ)對θ向量求一階導數，使得J』(θ)=0的時候J(θ)便取到了最值，同時也就可以求出相應的θ向量了，最終通過求解可以得出

其中X為數據矩陣，由於今天是針對logistic回歸進行求解，線性回歸的部分我會一筆帶過，有興趣的同學可以看http://blog.csdn.net/fishmemory/article/details/51586203這個博客，其中對於線性回歸的求導講解的比較詳細。

雖然用上述的方法可以直接求出θ向量，但一般我們並不會用這種方式去求解，因為矩陣的逆運算非常消耗計算機的資源，吳恩達建議特徵維度小於10,000時可以用直接求解的方式，但是大於10,000時要改用Gradient Descent方法，那麼第二個問題來了，什麼是Gradient Descent？Gradient Descent也就是常說的梯度下降，我們再來舉一個小例子：

對於 $y=x^{2}$ 這個函數，它開口向上存在最小值，取最小值時x=0，這個非常簡單，但如果我的起始點選在了A點，我應該怎樣讓A移動到x=0，取到y的最小值呢？這時候就要用到梯度下降的方法，我們知道如果用y對x求導，並且取 $x = x_{A}$ 可以求出A這個點的斜率，斜率是一個有方向的量，如果我們取 $x_{A} = x_{A} -varepsilon imes(k_{A})$ （其中ε為學習效率是一個很小的數，一般可以取0.01， $k_{A}$ 為y在 $x = x_{A}$ 時的導數值），如果一直不斷的迭代更新最終A就會慢慢的移動到x=0這一點，這也就完成了梯度下降，y也就取到了最小值，有同學肯定要問為什麼是要減去 $varepsilon imes(k_{A})$ ，以這個問題為例，A點的斜率是負數，A要走到最小值的取值位置要向斜率的反方向走才可以，因此要用減號。

一般的線性擬合問題的cost方程可以用梯度下降的方式求解權重的最優解，可以寫成：

其中α為學習效率，可以選擇一個很小的數。

至此我很粗略的講解了機器學習中的兩個很重要的概念，1、cost方程；2、Gradient Descent。一個是評判模型建立好壞的標準之一，另一個是求解模型中權重值θ的重要方法，有了這兩個基礎我們就可以來講解logistic回歸中的cost方程以及權重值θ是如何做Gradient Descent的了，那麼我們現在還是將目光投回到開篇講的logistic回歸中來吧。

Logistic回歸中cost方程以及梯度下降求解過程

Logistic回歸也存在著cost方程去評判模型建立的好壞，與一般線性回歸不同，他的cost方程如下（為了保證能與書中的結果匹配，這裡的表達式都與書中保持了一致）：

其中π為預測的結果：

如何去理解這個cost方程呢，這與線性回歸中 $sum_{a}^{b}{x}$ 這個符號表示累加不同， $prod_{a}^{b}$ 這個符號表示的是累乘，舉個例子如果現在就兩個客戶一好一壞，第一個如果是好客戶，那麼 $y_{1} = 0$ ，如果 $pi_{1} = 0.2$ 的話整個，對於第一個客戶他的cost就可以寫成 $0.2^{0} imes(1-0.2)^{1}$ ，第二個如果是壞客戶，那麼 $y_{2}=1$ ，如果 $pi_{2}=0.9$ 的話整個，對於第二個客戶他的cost就可以寫成 $0.9^{1} imes(1-0.9)^{0}$ ，如果有n個客戶將這些單個的cost結果進行累乘，不知道大家發現了什麼，其實這個cost方程就是對於每個獨立事件進行好壞的預測，並且這n個事件都預測對的概率，理論上來講如果模型是完美的那麼L=1，但這是做不到的，和線性問題一樣，我們要求cost的最小值，因此這個時候還需要將cost方程轉換一下：

$J(eta)=-logL$

對L做自然對數（如果沒有特殊說明log就代表ln），並求相反數，這樣就從原來求L的最大值，變成了求-logL的最小值了。接下去就與線性回歸一樣對θ求導做梯度下降便可，即：

其中β為模型的權重， $ar{I}^{-1}$ 為學習效率，g為J(β)對於β一階導數，也就是梯度的方向。

值得注意的是這時候《信用風險評分卡研究》的第七章中出現了一個小錯誤，這很有可能導致大家出現後續的迷惑，在第七章的115頁，書中定義g的演算法為：

然而到了書的116頁又變成了如下說法：

通過後續的計算我相信第二種說法才是正確的。

書中經過求導以後直接給出：

這對於第一次接觸logistic回歸的同學來說這實在是太難理解了，現在就讓我一步步對g進行求解吧，通過剛才的講解我們知道此時的cost方程已經變成： $J(eta)=-logL$ ，通過求解可以

得出：

其中：

這個時候我們可以對每一行數據進行求解，並將問題化簡為只有一個參數β，那這個問題就變得非常簡單：

其中 $L_{i}$ 表示某一個數據點的cost的方程，最終只要對結果進行累加即可，因此可以得出：

（這個地方不懂的同學可以去百度一下複合函數求導以及lnx的求導公式）

進一步化簡：

再次化簡：

最終整理一下：

是不是感覺和下邊的公式很像了：

沒錯，這個時候再考慮一下每一行數據還會有j個features並且整個數據集存在著i行數據最終便可以得到如下公式：

因此採用梯度下降的方法我們便可以求出最大似然估計分析表中的第一部分參數Estimate，也就是模型的權重參數。

讓我們再次回到這個公式：

其中的 $ar{I}^{-1}$ 是不是隨機選取的一個很小的值呢，在logistic回歸中並不是這樣，它等於-logL對β的二階導數，此時書中又犯了一個錯誤，在116頁：

讓我們再回到剛才的公式：

這個公式中，如果再一次對β求導可以得出（利用y=u/v求導：y=(uv-uv)/v^2這個公式）：

與剛才一樣我們將數據橫向以及縱向進行擴展最終得到：

進一步整理得到：

其中：

那麼這個信息矩陣和咱們的似然估計分析表有什麼關係呢？當然有，因為其中的Standard Error就是根據 $ar{I}^{-1}$ 求得的，書中提出，Standard Error就是 $ar{I}^{-1}$ 的對角元素的平方根，什麼意思？我們還是通過一段python代碼來演示一下把（代碼以及數據集我會放在github中供大家練習）：

最大似然估計參數求解

第一步導入一批數據：

第二步利用python中statsmodels模塊求logistic回歸：

觀察似然估計結果：

此時我們將會得到一個和sas當中十分相似的結果，其中第一列的coef實際上就是sas中求出的Estimate參數。

第三步，根據公式我們知道：

因此在python中輸入如下代碼：

最終輸出結果如下：

對比剛才的輸出結果，的確保持了一致：

好的，這個時候我們已經解決了最大似然估計表中的兩個參數了，接下來我們來看看Chi-Square這個值怎麼算出來的，大家可以看出python出現的z值這個統計量，那其中有沒有Chi-Square這個統計量呢，答案是肯定的。

第四步，求出卡方統計量：

大家可以看出實際上卡方值就是z值的平方，比如 $(-39.426)^2≈1554.409$ ，那麼這個卡方統計量又是怎麼得到的呢，很簡單這個值等於 $ChiSquare=(coef/std err)^2$

例如intercept的卡方值等於 $(-3.5929/0.091)^2 ≈ 1558.86$ ，大家可以驗證一下是否所有的卡方值都滿足這個公式：

接下去我們再來看看Pr>ChiSq這個值是怎麼求的吧，p值檢驗是logistic回歸最後輸出結果中一個十分重要的統計量，一般要求變數的p值不大於0.05，那麼這個p值是怎麼求解的呢，其實這個p值不是求出來的，而是查表查出來的：

通過求解我們發現這些參數的卡方值都是服從自由度為1的卡方分布的，因此對比書中的TmAtAddress這個變數的p值我們可以看出當卡方值約為0.0095時p值大約在0.9-0.95之間，其他的值大家也可以一一驗證一下：

此時我們也可以得出一個推論，那就是只要卡方值大於3.84，那麼這個變數基本上滿足了p值的要求。

最後一個問題，這個95%的置信區間又是怎麼求來的呢？我們知道沃爾德置信區間也被稱為正態分布區間，因此logistic回歸求解出的參數滿足如下公式：

其中 $sigma$ 就是求解出的std err值，n=1，通過查表得出 $mu_{0.025}$ =1.96，通過python代碼驗證結果如下：

對比統計結果可以看出該結果完全吻合：

同樣我們來驗證一下sas中的結果：

以Intercept為例：1.56+0.3598*1.96=2.265，1.56-0.3598*1.96=0.8547，結果基本上也吻合。

《信用風險評分卡研究》中最大似然估計分析表中的參數如何求解都已經講完了，這些內容都是我這段時間自己總結的一些小經驗，因此過程中難免會有一些講錯或者理解錯的地方，希望看到的同學能夠指正，尤其是其中學習效率為什麼選擇信息矩陣以及std err求解的原理我還沒有想清楚，希望知道原理的同學能給大家分享一下。

最後測試的代碼和數據集可以在這個鏈接中下載：http://pan.baidu.com/s/1dE3MyWD，希望大家能在本次的講解中有所收穫！