《信用风险评分卡研究》中最大似然估计分析表的解读

《信用风险评分卡研究》中最大似然估计分析表的解读

众所周知,Logistic回归在信用评分卡开发过程中起到了核心的作用,其回归的结果可以转换成一个汇总表,形成所谓的标准评分卡,除此之外在《信用风险评分卡研究》一书中还提出了大量对于模型评判的指标,验证模型的效果,可以说这些指标不仅保证了模型的解释性也为模型提供了良好的统计学基础,今天我想谈一谈的就是模型评判的一些重要指标—最大似然估计分析表,我们先看一下书中第七章出现的两张表:

  • 了解cost方程&梯度下降

不知道大家有没有和我有着一样的疑惑,那就是这些Estimate,Standard Error,Chi-Square,Pr>ChiSq这些参数到底是什么以及这些参数的怎么算出来的,他们的含义又都是什么。今天我就将从一个初学者的角度来解读一下logistic回归中这些指标到底怎么算出来的。

首先我们还是先回顾一下logistic回归的公式:

其中 \beta_{i} 代表每一个x对应的权重,而这个公式最终就是要通过这些x变量和对应的权重去预测最终一行数据输入后其等于1的概率,通过化简最终得出:

其中z = \beta^{T} x,β和x都是向量,那么如果当我们最终求出了一个P值以后我们怎么评判这个结果的好坏呢,这个时候我们一般会使用cost方程,好了第一个问题,什么是cost方程,让我们先跳出这个logistic回归,回到最简单的线性回归问题,我们先看一个实际中的简单例子:

对于如图所示的一堆散点,最终我们得到了一条红色的拟合曲线 h_{\theta} = \theta_{0}+\theta_{1} x,我们一般会使用如下方式去评判曲线拟合的好坏:

其中J(θ)为拟合的cost方程, h_{0}(x^{(i)}) 表示预测值, y^{(i)} 表示实际值,如果拟合的越好那么J(θ)的值就会越小,我们做机器学习的目的就是要找出合适的 \theta_{0}\theta_{1} 使得这个J(θ)最小,学习过导数的同学都知道,只要让J(θ)对θ向量求一阶导数,使得J’(θ)=0的时候J(θ)便取到了最值,同时也就可以求出相应的θ向量了,最终通过求解可以得出

其中X为数据矩阵,由于今天是针对logistic回归进行求解,线性回归的部分我会一笔带过,有兴趣的同学可以看blog.csdn.net/fishmemor这个博客,其中对于线性回归的求导讲解的比较详细。

虽然用上述的方法可以直接求出θ向量,但一般我们并不会用这种方式去求解,因为矩阵的逆运算非常消耗计算机的资源,吴恩达建议特征维度小于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\times(k_{A}) (其中ε为学习效率是一个很小的数,一般可以取0.01, k_{A} 为y在 x = x_{A} 时的导数值),如果一直不断的迭代更新最终A就会慢慢的移动到x=0这一点,这也就完成了梯度下降,y也就取到了最小值,有同学肯定要问为什么是要减去 \varepsilon\times(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}\times(1-0.2)^{1} ,第二个如果是坏客户,那么 y_{2}=1,如果 \pi_{2}=0.9 的话整个,对于第二个客户他的cost就可以写成 0.9^{1}\times(1-0.9)^{0} ,如果有n个客户将这些单个的cost结果进行累乘,不知道大家发现了什么,其实这个cost方程就是对于每个独立事件进行好坏的预测,并且这n个事件都预测对的概率,理论上来讲如果模型是完美的那么L=1,但这是做不到的,和线性问题一样,我们要求cost的最小值,因此这个时候还需要将cost方程转换一下:

J(\beta)=-logL

对L做自然对数(如果没有特殊说明log就代表ln),并求相反数,这样就从原来求L的最大值,变成了求-logL的最小值了。接下去就与线性回归一样对θ求导做梯度下降便可,即:

其中β为模型的权重, \bar{I}^{-1} 为学习效率,g为J(β)对于β一阶导数,也就是梯度的方向。

值得注意的是这时候《信用风险评分卡研究》的第七章中出现了一个小错误,这很有可能导致大家出现后续的迷惑,在第七章的115页,书中定义g的算法为:

然而到了书的116页又变成了如下说法:

通过后续的计算我相信第二种说法才是正确的。

书中经过求导以后直接给出:

这对于第一次接触logistic回归的同学来说这实在是太难理解了,现在就让我一步步对g进行求解吧,通过刚才的讲解我们知道此时的cost方程已经变成: J(\beta)=-logL ,通过求解可以

得出:

其中:

这个时候我们可以对每一行数据进行求解,并将问题化简为只有一个参数β,那这个问题就变得非常简单:

其中 L_{i} 表示某一个数据点的cost的方程,最终只要对结果进行累加即可,因此可以得出:

(这个地方不懂的同学可以去百度一下复合函数求导以及lnx的求导公式)

进一步化简:

再次化简:

最终整理一下:

是不是感觉和下边的公式很像了:


没错,这个时候再考虑一下每一行数据还会有j个features并且整个数据集存在着i行数据最终便可以得到如下公式:

因此采用梯度下降的方法我们便可以求出最大似然估计分析表中的第一部分参数Estimate,也就是模型的权重参数。

让我们再次回到这个公式:

其中的 \bar{I}^{-1} 是不是随机选取的一个很小的值呢,在logistic回归中并不是这样,它等于-logL对β的二阶导数,此时书中又犯了一个错误,在116页:

让我们再回到刚才的公式:

这个公式中,如果再一次对β求导可以得出(利用y=u/v求导:y'=(u'v-uv')/v^2这个公式):

与刚才一样我们将数据横向以及纵向进行扩展最终得到:

进一步整理得到:

其中:

那么这个信息矩阵和咱们的似然估计分析表有什么关系呢?当然有,因为其中的Standard Error就是根据 \bar{I}^{-1} 求得的,书中提出,Standard Error就是 \bar{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求解的原理我还没有想清楚,希望知道原理的同学能给大家分享一下。

最后测试的代码和数据集可以在这个链接中下载:pan.baidu.com/s/1dE3MyW,希望大家能在本次的讲解中有所收获!

编辑于 2017-10-30 23:03