python3 逻辑回归（含推导）

最新推荐文章于 2022-04-14 10:44:15 发布

我不爱机器学习

最新推荐文章于 2022-04-14 10:44:15 发布

阅读量1.2k

点赞数 3

分类专栏： python 数据分析

本文链接：https://blog.csdn.net/mengjizhiyou/article/details/103117274

版权

数据分析同时被 2 个专栏收录

25 篇文章 12 订阅

订阅专栏

python

24 篇文章 5 订阅

订阅专栏

1、比率（odds）、概率（probability，可能性）和似然（likelihood)

例如：成功7次，失败5次

比率 =某一类样本数/另一类样本数}=7/5
概率 = 某一类样本数/所有可能情况的样本=7/（7+5）
如果用theta表示环境参数，y表示结果，则概率(可能性)表示为： $P(y|\theta )$ (条件概率,已知theta），theta称为前置条件；在已知y的条件下求theta则称为似然： $P(y|\theta)$ ，一般采用极大似然法求参数theta，再用梯度下降法和牛顿法等求其最优解。

注：梯度下降法使用条件是该函数是凸函数，单调可微，连续

2、伯努利分布与二项分布

伯努利分布（Bernoulli distribution)又叫两点分布或0-1分布。记成功概率为p,失败概率为q=1-p.其中0<p<1.公式：

也可以写成 $f(y|p)=p^y(1-p)^{1-y}$ $y\in {(0,1)}$

二项分布是n 个独立的伯努利分布，连乘。

3、极大似然法（maximum likelihood method)

主要用于求解参数
多个独立事件的似然函数： $L=\prod_{i=1}^Np_i$ (N为样本数）
对多项乘积的求导往往非常复杂，但是对于多项求和的求导却要简单的多，对数函数不改变原函数的单调性和极值位置(对数化似然函数（log-likelihood），而且根据对数函数的性质可以将乘积转换为加减式，这可以大大简化求导的过程： $log(L)=log(\prod_{i=1}^Np_i)=\sum _{i=1}^Nlog(p_i)$

4、逻辑回归

用于分类，不是回归

4.1 联系函数（link function)

线性回归表达式： y = w^Tx+b ；

逻辑回归对y求对数（对数线性转化：log-linear regression)： ln(y) = w^Tx+b or y =exp( w^Tx+b) ,实际上是试图让接近y；

一般对于任何单调函数g(·）可将非线性函数转化为线性函数并且单调可微，令：

g(y) = w^Tx+b or $y = g^{-1}(w^Tx+b)$ 则该模型称为广义线性模型（generalized linear model),g(·)称为联系函数(link funjction).对数线性回归是广义模型在g(·）=In(·）时的特例。

4.2 单位阶跃函数（unit-step function）

考虑二分类任务，输出标记 $y\in (0,1)$ ,，而线性回归模型产生的预测值 ${\color{Red} z = w^Tx+b}$ 是实值。于是，将z转化为0/1值。最理想是单位阶跃函数:

$y=\left\{\begin{matrix} 0, & z<0\\ 0.5, &z=0 \\ 1,& z>0 \end{matrix}\right.$

即预测值z大于零就判为正例，小于零判为反例，预测值的临界值则可任意判别，如下图所示。

单位阶跃函数是不连续的，因此不可以用作 $g^{-1}(\cdot )$ ,因此需要找到能在一定程度上近似单位阶跃函数的替代函数（surrogate function),并且单调可微。对数几率函数（logistic function)正是这样一个函数： ${\color{Magenta} y = \frac{1}{1+e^{-z}}}$

对数几率函数是一种‘Sigmoid函数’,S型

4.3 逻辑回归推导

根据对数几率函数可得：

$y = \frac{1}{1+e^{-(w^Tx+b)}}\Rightarrow ln\frac{y}{1-y}=w^Tx+b$

上式中，若将y视为样本x作为正例的可能性，则1-y是其反例可能性，两者的比值= $\frac{y}{1-y}$ ，称为几率（odds)，反映了x作为正例的相对可能性。对几率取对数则得到对数几率（log odds,亦称logit)。

将y视为类后验概率估计 p(y=1|x) , p(y=0|x) ,则上式可写为：

$ln\frac{p(y=1|x)}{p(y=0|x)}=w^Tx+b$

$ln\frac{p(y=1|x)}{1-p(y=1|x)}=w^Tx+b\Rightarrow p(y=1|x)=\frac{e^{w^Tx+b}}{1+e^{w^Tx+b}}$

$p(y=0|x)=1-p(y=1|x)\Rightarrow p(y=0|x)=\frac{1}{1+e^{w^Tx+b}}$

于是，可通过极大似然法估计w,b，对率回归模型最大化对数似然(以e为底数），i代表第i个样本：

$L(w,b)=\prod_{i=1}^Np(y_i=1|x)^{y_i}p(y_i=0|x)^{1-yi}$

$l(w,b)=log(L(w,b))=log(\prod_{i=1}^Np(y_i=1|x)^{y_i}p(y_i=0|x)^{1-yi})$

$=\sum _{i=1}^Nlog(p(y_i=1|x)^{y_i}p(y_i=0|x)^{1-yi}))$

$=\sum _{i=1}^N[log(p(y_i=1|x)^{y_i})+log(p(y_i=0|x)^{1-yi}))]$

$=\sum _{i=1}^N[y_ilog( \frac{e^{w^Tx_i+b}}{1+e^{w^Tx_i+b}})+(1-y_i)log(\frac{1}{1+e^{w^Tx_i+b}})]$

$=\sum _{i=1}^N[y_ilog(e^{w^Tx_i+b})-y_ilog(1+e^{w^Tx_i+b})+(y_i-1)log(1+e^{w^Tx_i+b})]$

$=\sum _{i=1}^N[y_ilog(e^{w^Tx_i+b})-y_ilog(1+e^{w^Tx_i+b})+y_ilog(1+e^{w^Tx_i+b})-log(1+e^{w^Tx_i+b})]$

$=\sum _{i=1}^N[y_ilog(e^{w^Tx_i+b})-log(1+e^{w^Tx_i+b})]$

$=\sum _{i=1}^N[y_i(w^Tx_i+b)-log(1+e^{w^Tx_i+b})]$

对 l(w,b) 求极大值，得到w,b的估计。

为了方便，将w,x扩充为 w=(w_1,w_2,...,w_n,b)^T , x=(x_1,x_2,...,x_n,1)^T ；在求参数时， l(w,b) 最大化，就是对 -l(w,b) 最小化，因此 l(w)=-l(w,b) (加负号是为了从对数似然函数最大化问题等价变化到代价函数最小化问题)

求w的一次导数： $\bigtriangledown_w l(w)=\frac{\partial l(w)}{\partial w}=-\frac{\partial \sum _{i=1}^N[y_i(wx_i)-log(1+e^{wx_i})]}{\partial w}$

$=-\sum _{i=1}^N(y_ix_i-\frac{x_ie^{wx_i}}{1+e^{wx_i}})$

$=-\sum _{i=1}^Nx_i(y_i-\frac{e^{wx_i}}{1+e^{wx_i}})$

$=-\sum _{i=1}^Nx_i(y_i-g(wx_i))$

or $=-\sum _{i=1}^Nx_i(y_i-p(y=1|x_i))=-X^T(Y-\widehat{Y})$

根据梯度定义写成矩阵形式：

$\bigtriangledown_w l(w)=\begin{bmatrix} \frac{\partial l(w)}{\partial w_0}\\ \frac{\partial l(w)}{\partial w_1} \end{bmatrix} =\begin{bmatrix} \sum _{i=1}^Nx_0(g(wx_i)-y_i)\\ \sum _{i=1}^Nx_1(g(wx_i)-y_i) \end{bmatrix} =\begin{bmatrix} x_0^T(g(wx)-y)\\ x_1^T(g(wx)-y) \end{bmatrix}$

=X^T(g(WX)-Y)

令 $\bigtriangledown l(w)$ =0，发现没有解析解，只能用数值法（梯度下降，牛顿迭代...）

求w的二次导数：若采用梯度下降则不需计算二阶导，牛顿迭代需要计算二阶导

$\bigtriangledown_{ww} l(w)=\frac{\partial^2 l(w)}{\partial w\partial w^T}=\frac{\partial-\sum _{i=1}^Nx_i(y_i-\frac{e^{wx_i}}{1+e^{wx_i}})}{\partial w\partial w^T}$

$=\frac{\partial \sum _{i=1}^N\frac{e^{wx_i}}{1+e^{wx_i}}}{\partial w\partial w^T}=\frac{\partial \sum _{i=1}^Nx_i\frac{e^{wx_i}+1-1}{1+e^{wx_i}}}{\partial w\partial w^T}=\frac{\partial \sum _{i=1}^Nx_i(1-\frac{1}{1+e^{wx_i}})}{\partial w\partial w^T}$

$=\frac{\partial \sum _{i=1}^N-\frac{x_i}{1+e^{wx_i}}}{\partial w\partial w^T}$

$=\sum_{i=1}^N( -1*(-x_i)*(1+e^{wx_i})^{-2}*e^{wx_i}*x_i)$

$=\sum_{i=1}^Nx_i^2\frac{e^{wx_i}}{(1+e^{wx_i})^{2}}=\sum_{i=1}^Nx_i^2\frac{e^{wx_i}\cdot 1}{(1+e^{wx_i})^{2}}$

$=\sum_{i=1}^Nx_i^2\frac{e^{wx_i}\cdot 1}{(1+e^{wx_i})(1+e^{wx_i})} =\sum_{i=1}^Nx_i^2\frac{e^{wx_i}}{(1+e^{wx_i})}\cdot \frac{1}{(1+e^{wx_i})}$

$=\sum_{i=1}^Nx_ix_i^Tp(y=1|x_i)p(y=0|x_i)=X^TP(1-P)X=X^TWX$

数值求解可以用于参数求解（可令偏导数=0，但矩阵方法不可解，所以采用数值求解）：

采用牛顿迭代法： $\frac{y-f(x_0)}{x-x_0}={f}'(x_0)\Rightarrow y={f}'(x_0)(x-x_0)+f(x_0)$

$\bigtriangledown_{w} l(w^{t+1}) =\bigtriangledown_{ww} l(w)(w^{t+1}-w^{t})+\bigtriangledown_{w} l(w)=0$

$w^{t+1}=w^{t}-\frac{\bigtriangledown_{w} l(w)}{\bigtriangledown_{ww} l(w)}$ (Hessian Matrix)

将 $\bigtriangledown_{ww} l(w)=X^TWX$ ， $\bigtriangledown_{w} l(w)=-X^T(Y-\widehat{Y})$ ，带入Hessian Matrix:

$w^{t+1}=w^{t}+(X^TW^{(t)}X)^{-1}\cdot X(Y-\widehat{Y}^{(t)})$

最后将 $w^{t+1}$ 带入 ${\color{Magenta} p(X) = \frac{1}{1+e^{-wx}}}$ ,计算出概率，通过设定阈值进行分类，采取相关指标进行评价。

补充：

解析解与数值解的区别：

解析解(closed-form solution)：是指可以通过严格的公式推导出的精确解，如y=ax+b，只需要知道两个点就可以求出参数a，b。
数值解(numerical solution)：对于多个数据(n>>2)的情况，这时很有可能不存在解析解，我们只能借助数值方法去优化出一个近似的数值解。
之所以叫优化，是因为计算机的计算速度非常快，我们可以借助强大的计算能力去多次搜索和试错。从而一步步降低误差L。最简单的优化方法就是暴力搜索或随机试验，若找出最合适的a、b，则可以从(部分)实数空间中随机采样任意的a、b，并计算出对应模型的误差值L，然后从测试过的L中挑出最好的L，它所对应的a、b就可以作为要找的最优a、b。

4.4 逻辑回归正则化

l1 正则化：

l2 正则化：

弹性网络正则化：

4.5 多项逻辑回归

上面的逻辑回归模型是二项分类模型，适用于二分类，可将其推广为多项逻辑回归模型（multi-nominal logistic regression model)，用于多分类。假设离散型随机变量Y的取值集合是{1，2，...，K}，则多项逻辑回归模型是：

$\bg_white P(Y=k|x)=\frac{exp(w_kx)}{1+\sum_{k=1}^{k-1}exp(w_k\cdot x)}$ ， k=1,2,...,K-1

$\bg_white P(Y=K|x)=\frac{1}{1+\sum_{k=1}^{k-1}exp(w_k\cdot x)}$

4.6 python3

from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
import numpy as np
import matplotlib.pyplot as plt

X, y = load_iris(return_X_y=True)
X = X[:, :2]  # we only take the first two features
logreg = LogisticRegression(C=1e5, solver='lbfgs', multi_class='multinomial').fit(X, y)

# Plot the decision boundary. For that, we will assign a color to each
# point in the mesh [x_min, x_max]x[y_min, y_max]
x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
h = .02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

# put the result into a color plot
Z = Z.reshape(xx.shape)
plt.figure(1, figsize=(4, 3))
plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)

# plot also the training points
plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k', cmap=plt.cm.Paired)
plt.xlabel('Sepal length')
plt.ylabel('Sepal width')
plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xticks(())
plt.yticks(())
plt.show()