逻辑回归（logistic regression）

来源：互联网发布：tv365网络电视安卓版编辑：程序博客网时间：2024/05/16 14:06

一、数学上的逻辑回归

前面提到，逻辑回归是判别分析方法来分类的，即通过给定的数据x, 来直接得到其后验概率。且它得到的是线性分类边界。

回顾在贝叶斯准则中，利用0-1损失进行分类时，我们做法是以最大的后验概率的类 k，来作为依据。

${\hat G}(x) =\mathop {\arg \max }\limits_kPr(G = k|X=x)$

从而第k 类和第 l 类的分类边界通过使其概率相等来决定：即样本 x 在第k 类和第l 类有相等的后验概率。

$Pr(G = k|X=x) = Pr(G = l|X=x)$

如果我们我们对两个概率进行相除，并且取 log 的话，这样得到一个比率，在上述情况有

$log \frac{Pr(G = k|X=x)} {Pr(G = l|X=x)} = 0$

如果我们想要强制得到一个线性分类边界，我们可以假设这个函数可以线性表示，即

$log \frac{Pr(G = k|X=x)} {Pr(G = l|X=x)} = a_0^{(k,l)} + \sum \limits_{j=1}^p{a_j^{(k,l)}} x_j$

逻辑回归就是在这样一个假设的的基础上得到的。其中对不同的两类类，有不同的系数上标（k,l）, 其系数也不同。

这样看着很容易理解，那肯定也会有这样的疑问，那岂不是每两类就要找一个 $a_0^{(k,l)}$ ，那计算量岂不是很大！

然而在逻辑回归中，我们并不用每两类都要找一组系数，对于 K 类，我们只需要进行 K-1 对配对，找K-1 组系数就可以了。

假设开始

现在我们有 K 类，我们让第K 类（可以是任何一个类）作为一个基类，这样对于剩余的 K-1 类，得到 K-1 组系数情况：

$\begin{align*} log \frac{Pr(G = 1|X=x)} {Pr(G = K|X=x)}& = \beta_{10} + \beta_1^T x\\ log \frac{Pr(G = 2|X=x)} {Pr(G = K|X=x)}& = \beta_{20} + \beta_2^T x\\ & \vdots \\ log \frac{Pr(G = K-1|X=x)} {Pr(G = K|X=x)}& = \beta_{(K-1)0} + \beta_{K-1}^T x \end{align*}$

我们不必找到每对的分类边界，只要找到每个类和基类的分类边界的系数。

一旦找到了 K -1 类的 log 比率，对于任何一对（k, l）的log 比率，我们就可以通过上述得到的系数组推导出来，比如:

这样我们总共参数的个数为（K-1）×（p+1）其中(k-1) 是 k-1 组系数，每组中有p 维，再添加一个常数项 1。

为了方便，系数表示为: $\theta = \{ \beta_{10},\beta_1,\beta_{20},\beta_2,...,\beta_{(K-1))0},\beta_{K-1}\}$ , 其中后验概率的log 比率叫做 log-odds 或者logit transformations.

通过上述假设，我们得到以下两个后验概率公式：

$Pr(G=k|X=x)=\frac {exp(\beta_{k0}+\beta_k^T x)}{1+\sum \limits_{l=1}^{K-1}exp(\beta_{l0}+\beta_l^T x)} \\ \quad for \;k = 1,2,...,K-1$

$Pr(G=K|X=x)=\frac {1}{1+\sum \limits_{l=1}^{K-1}exp(\beta_{l0}+\beta_l^T x)}$

第一个公式，可以理解为，将Pr(G=K|X =x) 当作单位1，剩下的每个K-1 个与他比，其所占的部分是 $exp(\beta_{k0}+\beta_k^T x)$ ，因此对于一个 k 类，所占的比重为，自己的部分比上所有加和。

第二个公式也可以这样理解，还可以理解为概率的总和是 1，因此其对于K 类，其结果为第二个公式。

因此这样就得到了，众说的sigmoid 函数，具体得到就是这样了。（来自wiki）

逻辑回归比线性回归相比，它对于 x 是非线性函数，且概率在0-1之间，总和为 1。

二、参数估计

在得到上述参数后，要计算概率，就要求解参数 $\theta$ ，然后最大化条件似然G，在给定x 之后，我们感兴趣的不是 X 的分布情况，而是在哪类的概率最大。给一点 xi, 其属于 k 类的后验概率表示为： $p_k(x_i;\theta)=Pr(G = k|X=x_i;\theta)$ 。对于第一个样本x1, 它属于g1 类，则其后验概率为

$Pr(G = g_1|X=x_1)$

因为样本是独立的，因此，这些N 个样本点每个类为 gi, 所以其后验概率为。即联合条件似然就是其条件概率的乘积。

$\prod\limits_{i=1}^{N} {Pr(G = g_i|X=x_i)}$

对参数的似然估计即为 log条件似然，通过使其最大化得到。

$\begin{align*} l(\theta) &= \sum \limits_{i=1}{N}{logPr(G = g_i|X=x_i)}\\ &= \sum \limits_{i=1}{N}{logp_{g_i}(x_i;\theta)} \end{align*}$ 、

具体计算过程参见下面二分类实战部分。

三、二分类实战

先看二分类，对于多类是相似的情况。

对于二分类，如果 gi = class 1, 表示为yi = 1; 如果 gi = class 2, 表示为 yi = 0. 将类标签变为0 和1 使其更加简化。

对于上式log 似然中的概率项，用

$p_1(x;\theta)=p(x;\theta)$

$p_2(x;\theta)=1- p_1(x;\theta)=1- p(x;\theta)$

,因为这两者后验概率和为1。

因为 K = 2, 我们只有一个线性多项式，一个分类决策边界，参数列表为 $\theta = \{\beta_{10}, \beta_1\}$ , 现在我们在这里表示为 $\beta = $\beta_{10}, \beta_1$^T$ , 这是一个列向量。

现在我们要做的就是一步步，简化，计算，求解

简化

对于 yi = 1, 有 gi = 1 这时其 log 似然为

$\begin{align*} logp_{g_i}(x;\beta)&=logp_1(x;\beta) \\ &= 1\cdot logp(x;\beta)\\ &=y_i logp(x;\beta) \end{align*}$

同样对于 yi = 0, 有 gi = 2 , 其形式为：

$\begin{align*} logp_{g_i}(x;\beta)&=logp_2(x;\beta) \\ &= 1\cdot log(1-p(x;\beta))\\ &=(1-y_i) log(1-p(x;\beta)) \end{align*}$

因为 yi = 0 或者 1-yi = 0, 因此我们可以将以上两部分加在一起，得到

$logp_{g_i}(x;\beta)= y_i logp(x;\beta)+(1-y_i) log(1-p(x;\beta))$

这样做不用单独拿出来分析，对于不同的yi 都适用。相加的log 似然形式也可得到哦简化：

$l(\beta) = \sum \limits_{i = 1}^{N} logp_{g_i}(x_i;\beta) = \sum \limits_{i = 1}^{N}[y_ilogp(x_i;\beta)+(1-y_i)log(1-p(x_i;\beta))]$

在 $\beta = $\beta_{10}, \beta_1$^T$ 中有p+1 个参数

则参数和函数对应关系即 :

$\beta = \left( {\begin{array}{*{20}{c}} \beta_{10}\\ {\begin{array}{*{20}{c}} \beta_{11}\\ \beta_{12}\\ \vdots \end{array}}\\ \beta_{1,p} \end{array}} \right) \Leftrightarrow x = \left( {\begin{array}{*{20}{c}} 1\\ {\begin{array}{*{20}{c}} x_{,1}\\ x_{,2}\\ \vdots \end{array}}\\ x{,p} \end{array}} \right)$

计算

有前面的假设，哦们得到logistic 回归模型为：

$\begin{align*} p(x;\beta) &= Pr(G=1|X=x) = \frac {exp(\beta^Tx)}{1+exp(\beta^Tx)} \\ 1-p(x;\beta) &= Pr(G=2|X=x) = \frac {1}{1+exp(\beta^Tx)} \end{align*}$

则 log-似然结果：

$l(\beta) = \sum \limits_{i =1}^N[y_i\beta^Tx_i-log(1+e^{\beta^Tx_i})]$

似然函数求参数，即最大化似然函数，这样对 $\beta_{1j}$ 求导j =0,1,2...,p ，就好：

$\begin{align*} \frac {\partial l(\beta)}{\partial \beta_{1j}} &= \sum \limits_{i=1}^N{y_ix_{ij}-\sum \limits_{i=1}^N \frac {x_{ij}e^{\beta^Tx_i}}{1+e^{\beta^Tx_i}}} \\ &= \sum \limits_{i=1}^N{y_ix_{ij}} -\sum \limits_{i=1}^Np(x;\beta)x_{ij} \\ &= \sum \limits_{i=1}^Nx_{ij}({y_i}-p(x_i;\beta)) \end{align*}$

对上述式子进行矩阵运算即可，其中 xi 是列向量 yi 是标量。

为了解 p+1 个非线性等式 $\frac {\partial l(\beta)}{\partial \beta_{1j}} = 0, \; j = 0,1..p$ ，我们利用 Newton迭代法进行求解：

牛顿迭代法需要二阶导数 Hessian 矩阵：

$\frac {\partial^2 l(\beta)}{\partial \beta\partial \beta^T} = -\sum \limits_{i=1}^N {x_ix_i^T p(x_i;\beta)(1-p(x_i;\beta))}$

在第 j 行第 n 列的情况如下：

$\frac {\partial l(\beta)}{\partial \beta_{1j}\partial \beta_{1n}} \\ &= - \sum \limits_{i=1}^N \frac {(1+e^{\beta^Tx_i})e^{\beta^Tx_i}x_{ij}x_{in} - (e^{\beta^Tx_i})^2{x_{ij}x_{in}}} {(1+e^{\beta^Tx_i})^2}\\ &= - \sum \limits_{i=1}^N x_{ij}x_{in} p(x_i;\beta) - x_{ij}x_{in} p(x_i;\beta)^2 \\ &= -\sum \limits_{i=1}^N {x_ix_i^T p(x_i;\beta)(1-p(x_i;\beta))}$

这样我们利用牛顿迭代，不断更新参数：

$\beta^{new} = \beta^{old}-(\frac {\partial^2 l(\beta)}{\partial \beta\partial \beta^T})^{-1} \frac {\partial l(\beta)}{\partial \beta}$

这样，如果给了一个旧的参数 $\beta^{old}$ ，通过得到Hessian 矩阵与一阶导数的乘积，这样就可对 $\beta^{old}$ ，进行更新得到 $\beta^{new}$ .

求解

迭代的过程可用以下矩阵形式来进行：

y: 为 yi 的列向量来表示
X: 为输入矩阵 N×(p+1) 表示
p: 拟合概率的N维向量，第 i 个数据样本表示为 $p(x_i;\beta^{old})$
W: N×N为对角权重矩阵，第 i 个对角元素为 $p(x_i;\beta^{old})(1-p(x_i;\beta^{old}))$
在上述基础上构建一阶导数表示为：

$\frac {\partial l(\beta)}{\partial \beta} = X^T(y-p)$

Hessian 矩阵表示为：

$\frac {\partial^2 l(\beta)}{\partial \beta \partial \beta^T} = -X^TWX$

因而牛顿步更新之后为：

$\begin{align*} \beta^{new} &= \beta^{old}+ (X^TWX)^{-1}X^T(y-p) \\ &= (X^TWX)^{-1}X^TW(X\beta^{old}+W^{-1} (y-p))\\ &=(X^TWX)^{-1}X^TWz \end{align*} \\ where ~~ z = X\beta^{old}+W^{-1} (y-p)$

如果 z 看作是 X 输入下的一个响应输出，则 $\beta^{new}$ 求解过程是一个带权的最小二乘问题：

$\beta^{new} \leftarrow \mathop {\arg \min}\limits_\beta {(z-X\beta)^TW(z-X\beta)}$

回想想线性回归中的最小二乘是去解：

$\mathop {\arg \min}\limits_\beta {(z-X\beta)^T(z-X\beta)}$

z 做为一个适应性的响应，这就是 IRLS（iteratively reweighted least squares）算法。

其伪代码算法过程为：

因为在这里要求 W 的逆进行逆运算，因此计算效率有所下降，这里引入了修改后的高效运算方法，两者是等价的。

这里在计算时，要注意运算第三步，对每个求概率是都是求p ，没有求 1-p。

c++ 代码算法实现：

利用了男性、女性由身高体重进行区分的分类，总共210 个样本点，且最终训练数据误差率为11.9048% （后续可采用交叉验证来试验）。

/*@Newton-Method : @iteratively reweighted least squares@x_train : train sample of data@y_train : train result of data@weight : the train weight result.@use the Computational Efficiency@use the matrix term to express the process*/void IterReweightedLS(Matrix &x_train, Matrix &y_train, Matrix &weight){//change the y as 1->1, 2->0for (int row = 0; row < y_train.GetMatrixRow(); row++){if (y_train.GetMatrixValue(row, 0) == 2.0)y_train.SetMatrixOnePosition(row, 0, 0.0);}//use the pro-resadual to judge weather convergeMatrix odl_pro_use_for_converge(y_train.GetMatrixRow(), y_train.GetMatrixCol());odl_pro_use_for_converge.SetMatrixAllOne();odl_pro_use_for_converge = odl_pro_use_for_converge.MatrixMultiNumPerRow(-1, 0);int iter_count = IterNumber; while (iter_count--){// use the matrix to express the culculate// the probability of the sigmoid express formMatrix temp(x_train.GetMatrixRow(), 1);Matrix probability(y_train.GetMatrixRow(), y_train.GetMatrixCol());temp = x_train.MatrixMulti(weight);//another thing the probability all be exp(xw)/exp(xw)+1, because we have connected it all togetherfor (int pro_row = 0; pro_row < y_train.GetMatrixRow(); pro_row++){double pro = sigmoid(temp.GetMatrixValue(pro_row, 0));probability.SetMatrixOnePosition(pro_row, 0, pro);}//the W of diagonal matrix use the value p*(1-p)Matrix W(y_train.GetMatrixRow(), y_train.GetMatrixRow()); W.SetMatrixAllZero();Matrix one_cut_pro(probability.GetMatrixRow(),1);one_cut_pro.SetMatrixAllOne();Matrix pro_multi_one_cut_pro(probability);pro_multi_one_cut_pro = probability.MatrixPointMlutiMatrix((one_cut_pro.MatrixCut(probability)));W.SetDialogUsVector(pro_multi_one_cut_pro);//x_hat for culmulate efficient// remember the iter pre and after.Matrix x_hat = x_train;for (int row = 0; row < x_train.GetMatrixRow(); row++){x_hat = x_hat.MatrixMultiNumPerRow(pro_multi_one_cut_pro.GetMatrixValue(row, 0), row);}//culmulate the (x^T x_hat)^(-1)*x^T*(Y-P) // first LUP factorization// and then culmulate the parameter w_redisualMatrix x_trans = x_train.MatrixTrans();Matrix x_factorization = x_trans.MatrixMulti(x_hat);Matrix right_value = x_trans.MatrixMulti(y_train.MatrixCut(probability));//(x^T x_hat)^(-1) if the x_factorization//use the lup factorization to culmulateLUPFactorization matrix_factor(x_factorization);matrix_factor.Factorization();right_value = (matrix_factor.GetP()).MatrixMulti(right_value);//use the l-down triangle and u-up triange to compute itMatrix intermediray_vector = right_value.MatrixDivideDownTriangle(matrix_factor.GetL());Matrix terminal_vector = intermediray_vector.MatrixDivideUpTriangle(matrix_factor.GetU());weight = weight.MatrixPlus(terminal_vector);Matrix pro_redasual = probability.MatrixCut(odl_pro_use_for_converge);//whether converged  double sum_pro_resudal = pro_redasual.VectorL1Norm();//cout << sum_pro_resudal << endl;if (sum_pro_resudal < (NumberSamble * Epsilon)){cout << "the gradient l2 norm is: " << sum_pro_resudal << endl;cout << "the iter number is: " << IterNumber -  iter_count << endl;weight.ShowMatrix();break;}odl_pro_use_for_converge = probability;}}

最终结果：收敛，迭代次数，最终更新参数，误差率。

并且利用 opencv 将结果进行了显示：

分类来看，效果还是不错的。

正则化（l2）逻辑回归

即用MAP代替 MLE, 并且计算是在其基础上修改目标函数，梯度和 Hessian 矩阵。

贝叶斯逻辑回归

逻辑回归我们并没有考虑先验，如果加入先验的清况，这时，我们希望逻辑回归是高斯分布的，且其先验形式为 $p(w)=N(w|0,V_0)$ , 这样其后验概率也是高斯的，即

$p(w|D)=N(w|\hat w,H^{-1})$

整个近似的过程是 Laplace 近似：因为我们想要后验是正太分布，但是逻辑回归的估计式子是没有相应的共轭先验与其相对应的，这样，我们近似过程如下。

假设其后验为

$p(\theta|D)=\frac {1}{Z}e^{-E(\theta)}$

其中 $E(\theta)$ 是能量函数，等于负log非正太log 后验， $E(\theta) = -logp(\theta,D)$ ，其中 z = p(D) 是正太常值。这样我们对 $E(\theta)$ 进行 Taylor 级数展开：

$E(\theta) \approx E(\theta^*) +(\theta-\theta^*)g+\frac {1}{2}(\theta-\theta^*)^TH(\theta-\theta^*)$

这样我们取在 $\theta^*$ 处取得最小能量值，且其 g 梯度为0 ，这样我们可得到正太的后验分布，其中：

$\hat p(\theta|D) & \approx \frac{1}{Z}e^{-E(\theta^*)} exp[-\frac{1}{2}(\theta-\theta^*)^TH(\theta-\theta^*)] = N(\theta|\theta^*, H^{-1})\\ Z = p(D) & \approx \int \hat p(\theta|D) d\theta = e^{-E(\theta^*)}(2\pi)^{D/2}|H|^{-\frac{1}{2}}$

这样就可得到想要的高斯分布。

另《Bayesian Data Analysis》介绍对于 weakly information of prior 可考虑用 strdent-分布，然后用EM 进行参数的估计。也可进行贝叶斯经验运用。

四、多类情况

多类的情况，是与二类的情况相类似，计算量增加了许多。参见 https://onlinecourses.science.psu.edu/stat557/node/56

五、总结

通过直接的判别分析，不利用先验就可进行分类，并且在直接分类上，如果对模型未知时往往到达不错的效果。并且可以利用 L2正则化进行不理想情况的修正；在贝叶斯模型中，加入先验时，可以利用近似方法来得到想要的模型。

六、参考资料

1. https://onlinecourses.science.psu.edu/stat557/node/52

2. 《Machine Learning A Probabilistic Perspective》第八章 logistic regression

3. 《Bayesian Data Analysis》

0 0