对数几率回归详细推导及python实现

菜鸟一枚，有错误欢迎指正，但不要喷。

在线性回归中，要求解的模型是
(1)y=wTx+b y=\mathbf{w^Tx}+b \tag{1} y=wTx+b(1)
实际上线性回归模型在稍加修改后也可以求解其他问题，例如考虑模型
(2)y=ewTx+b y=e^{\mathbf{w^Tx}+b} \tag{2} y=ewTx+b(2)
两边同时取对数后有
(3)lny=wTx+b lny=\mathbf{w^Tx}+b \tag{3} lny=wTx+b(3)
在形式上仍然是线性回归，但实际上已经是输入空间到输出空间的非线性映射。更一般的，考虑单调可微函数g(⋅)g(\cdotp)g(⋅)，令
(4)y=g−1(wTx+b) y=g^{-1}(\mathbf{w^Tx}+b) \tag{4} y=g−1(wTx+b)(4)
这样得到的模型称为广义线性模型，其中函数称为g(⋅)g(\cdotp)g(⋅)“联系函数”。

接下来就可以利用广义线性模型来解决分类问题，也就是对数几率回归。

在二分类任务中，输出标记y∈{0,1}y \in \{0,1\}y∈{0,1}，而线性回归模型z=wTx+bz=\mathbf{w^Tx}+bz=wTx+b产生的预测值是实值，于是需要将实值zzz转换为0/1值，最理想的函数是单位越阶函数(亦称Heaviside函数)
(5)y={0,z&lt;00.5,z=01,z&gt;0 y=\begin{cases} 0, &amp;z &lt; 0 \\ 0.5, &amp;z=0 \\ 1, &amp;z&gt;0 \end{cases} \tag{5} y=⎩⎪⎨⎪⎧​0,0.5,1,​z<0z=0z>0​(5)
即若预测值大于0就判为正例，小于0则判为反例，为0则任意判别。而单位越阶函数不连续，因此不能用作g−1(⋅)g^{-1}(\cdotp)g−1(⋅)，因此便有了替代函数
(6)y=11+e−z y=\frac{1}{1+e^{-z}} \tag{6} y=1+e−z1​(6)
这是一个Sigmoid函数，即形似S型的函数。

Sigmoid函数能够将z值转换为0和1之间的y值，其输出值在z=0附近变化很陡，能够很好的近似单位越阶函数。

将Sigmoid函数作为g−1(⋅)g^{-1}(\cdotp)g−1(⋅)代入式(4)得到
(7)y=11+e−(wTx+b) y=\frac{1}{1+e^{-(\mathbf{w^Tx}+b)}} \tag{7} y=1+e−(wTx+b)1​(7)
若y≥0.5y \geq 0.5y≥0.5则预测为正例，反之预测为反例。上式可以变化为
(8)lny1−y=wTx+b ln\frac{y}{1-y}=\mathbf{w^Tx}+b \tag{8} ln1−yy​=wTx+b(8)
若将yyy看作样本x\mathbf{x}x作为正例的可能性，则1−y1-y1−y是其反例的可能性。两者的比值
(9)y1−y \frac{y}{1-y} \tag{9} 1−yy​(9)
称为“几率"，反应了样本x\mathbf{x}x作为正例的相对可能性，取对数则得到"对数几率"
(10)lny1−y ln\frac{y}{1-y} \tag{10} ln1−yy​(10)
这就是对数几率回归名称的由来。对数几率回归在得到分类结果的同时还能得到近似概率预测。

接下来介绍如何确定w\mathbf{w}w和bbb。如果将式(8)中的yyy看作类后验概率p{y=1∣x}p\{y=1 \mid \mathbf{x}\}p{y=1∣x}，式(8)可重写为
(11)lnp(y=1∣x)p(y=0∣x)=wTx+b ln \frac{p(y=1 \mid \mathbf{x})}{p(y=0 \mid \mathbf{x})}=\mathbf{w^Tx}+b \tag{11} lnp(y=0∣x)p(y=1∣x)​=wTx+b(11)
显然有
(12)p(y=1∣x)=e(wTx+b)1+e(wTx+b) p(y=1 \mid \mathbf{x})=\frac{e^{(\mathbf{w^Tx}+b)}}{1+e^{(\mathbf{w^Tx}+b)}}\tag{12} p(y=1∣x)=1+e(wTx+b)e(wTx+b)​(12)

(13)p(y=0∣x)=11+e(wTx+b) p(y=0 \mid \mathbf{x})=\frac{1}{1+e^{(\mathbf{w^Tx}+b)}}\tag{13} p(y=0∣x)=1+e(wTx+b)1​(13)

于是可以通过极大似然估计来估计w\mathbf{w}w，bbb，给定数据集{(xi,yi)}i=1m\{(\mathbf{x_i}, y_i)\}^m_{i=1}{(xi​,yi​)}i=1m​，对率回归模型最大化对数似然
(14)l(w,b)=∑i=1mlnp(yi∣xi;w,b) l(\mathbf{w}, b)=\sum_{i=1}^{m}lnp(y_i \mid \mathbf{x_i};\mathbf{w},b)\tag{14} l(w,b)=i=1∑m​lnp(yi​∣xi​;w,b)(14)
为方便书写，令β=(w,b)T\pmb{\beta}=(\mathbf{w},b)^Tβ​β​​β=(w,b)T，x^=(x,1)T\hat{\mathbf{x}}=(\mathbf{x},1)^Tx^=(x,1)T，则wTx+b\mathbf{w^Tx}+bwTx+b可简写为βTx\pmb{\beta}^T\mathbf{x}β​β​​βTx，再令p1(x^;β)=p(y=1∣x^;β)p_1(\hat{\mathbf{x}};\pmb{\beta})=p(y=1 \mid \hat{\mathbf{x}};\pmb{\beta})p1​(x^;β​β​​β)=p(y=1∣x^;β​β​​β)，p0(x^;β)=p(y=0∣x^;β)=1−p1(x^;β)p_0(\hat{\mathbf{x}};\pmb{\beta})=p(y=0 \mid \hat{\mathbf{x}};\pmb{\beta})=1-p_1(\hat{\mathbf{x}};\pmb{\beta})p0​(x^;β​β​​β)=p(y=0∣x^;β​β​​β)=1−p1​(x^;β​β​​β)，接下来对似然函数进行化简，p(yi∣xi;w,b)p(y_i \mid \mathbf{x_i};\mathbf{w},b)p(yi​∣xi​;w,b)可以写作：
p(yi∣xi^;β)={p1(x^;β)yi=1p0(x^;β)yi=0 p(y_i \mid \hat{\mathbf{x_i}};\pmb{\beta})=\begin{cases} p_1(\hat{\mathbf{x}};\pmb{\beta}) &amp;y_i=1\\ p_0(\hat{\mathbf{x}};\pmb{\beta}) &amp;y_i=0 \end{cases} p(yi​∣xi​^​;β​β​​β)={p1​(x^;β​β​​β)p0​(x^;β​β​​β)​yi​=1yi​=0​
取对数
lnp(yi∣xi^;β)=yilnp1(x^;β)+(1−yi)lnp0(x^;β) lnp(y_i \mid \hat{\mathbf{x_i}};\pmb{\beta})=y_ilnp_1(\hat{\mathbf{x}};\pmb{\beta})+(1-y_i)lnp_0(\hat{\mathbf{x}};\pmb{\beta}) lnp(yi​∣xi​^​;β​β​​β)=yi​lnp1​(x^;β​β​​β)+(1−yi​)lnp0​(x^;β​β​​β)
取指数
p(yi∣xi^;β)=p1(x^;β)yi⋅p0(x^;β)1−yi p(y_i \mid \hat{\mathbf{x_i}};\pmb{\beta})=p_1(\hat{\mathbf{x}};\pmb{\beta})^{y_i}\cdotp p_0(\hat{\mathbf{x}};\pmb{\beta})^{1-y_i} p(yi​∣xi​^​;β​β​​β)=p1​(x^;β​β​​β)yi​⋅p0​(x^;β​β​​β)1−yi​
代入式(14)得
l(β)=∑i=1mln(p1(x^;β)yi⋅p0(x^;β)1−yi)=∑i=1m(yilnp1(x^;β)+(1−yi)lnp0(x^;β))=∑i=1m(yi(lnp1(x^;β)−lnp0(x^;β))+lnp0(x^;β))=∑i=1m(yilnp1(x^;β)p0(x^;β)+lnp0(x^;β))=∑i=1m(yiβTx^+ln11+e(wTx+b))=∑i=1m(yiβTx^−ln(1+e(wTx+b))) \begin{aligned} l(\pmb{\beta})&amp;=\sum_{i=1}^{m}ln(p_1(\hat{\mathbf{x}};\pmb{\beta})^{y_i}\cdotp p_0(\hat{\mathbf{x}};\pmb{\beta})^{1-y_i})\\ &amp;=\sum_{i=1}^{m}(y_ilnp_1(\hat{\mathbf{x}};\pmb{\beta})+(1-y_i)lnp_0(\hat{\mathbf{x}};\pmb{\beta}))\\ &amp;=\sum_{i=1}^{m}(y_i(lnp_1(\hat{\mathbf{x}};\pmb{\beta})-lnp_0(\hat{\mathbf{x}};\pmb{\beta}))+lnp_0(\hat{\mathbf{x}};\pmb{\beta}))\\ &amp;=\sum_{i=1}^{m}(y_iln\frac{p_1(\hat{\mathbf{x}};\pmb{\beta})}{p_0(\hat{\mathbf{x}};\pmb{\beta})}+lnp_0(\hat{\mathbf{x}};\pmb{\beta}))\\ &amp;=\sum_{i=1}^{m}(y_i\pmb{\beta}^T\hat{\mathbf{x}}+ln\frac{1}{1+e^{(\mathbf{w^Tx}+b)}})\\ &amp;=\sum_{i=1}^{m}(y_i\pmb{\beta}^T\hat{\mathbf{x}}-ln(1+e^{(\mathbf{w^Tx}+b)})) \end{aligned} l(β​β​​β)​=i=1∑m​ln(p1​(x^;β​β​​β)yi​⋅p0​(x^;β​β​​β)1−yi​)=i=1∑m​(yi​lnp1​(x^;β​β​​β)+(1−yi​)lnp0​(x^;β​β​​β))=i=1∑m​(yi​(lnp1​(x^;β​β​​β)−lnp0​(x^;β​β​​β))+lnp0​(x^;β​β​​β))=i=1∑m​(yi​lnp0​(x^;β​β​​β)p1​(x^;β​β​​β)​+lnp0​(x^;β​β​​β))=i=1∑m​(yi​β​β​​βTx^+ln1+e(wTx+b)1​)=i=1∑m​(yi​β​β​​βTx^−ln(1+e(wTx+b)))​
于是最大化的目标为
l(β)=∑i=1m(yiβTx^−ln(1+e(wTx+b))) l(\pmb{\beta})=\sum_{i=1}^{m}(y_i\pmb{\beta}^T\hat{\mathbf{x}}-ln(1+e^{(\mathbf{w^Tx}+b)})) l(β​β​​β)=i=1∑m​(yi​β​β​​βTx^−ln(1+e(wTx+b)))
等价于最小化
(15)l(β)=∑i=1m(−yiβTx^+ln(1+e(wTx+b))) l(\pmb{\beta})=\sum_{i=1}^{m}(-y_i\pmb{\beta}^T\hat{\mathbf{x}}+ln(1+e^{(\mathbf{w^Tx}+b)}))\tag{15} l(β​β​​β)=i=1∑m​(−yi​β​β​​βTx^+ln(1+e(wTx+b)))(15)
这是关于β\pmb{\beta}β​β​​β的高阶可导凸函数，使用梯度下降法即可求解，仿照scikit-learn的实现模式，代码如下：

import numpy as np

# 对数几率回归
class LogitRegression:
def __init__(self):
self._beta = None
self.coef_ = None
self.interception_ = None

def _sigmoid(self, x):
return 1. / (1. + np.exp(-x))

def fit(self, X, y, eta=0.01, n_iters=1e4, eps=1e-8):
assert(X.shape[0] == y.shape[0])

def L(beta, X_b, y):
try:
vec = beta.dot(X_b.T)
return np.sum(-vec * y + np.log(1 + np.exp(vec))) / len(X_b)
except Exception:
return float('inf')

def dL(beta, X_b, y):
try:
vec = beta.dot(X_b.T)
coef =  np.exp(vec) / (1 + np.exp(vec))
return (-y.dot(X_b) + X_b.T.dot(coef)) / len(X_b)
except Exception:
return float('inf')

def gradient_decent(X_b, y, initial_beta, eta, n_iters, eps):
i_iters = 0
beta = initial_beta
while i_iters < n_iters:
gradient = dL(beta, X_b, y)
last_beta = beta
beta = beta - eta * gradient
if (abs(L(beta, X_b, y) - L(last_beta, X_b, y)) < eps):
break
i_iters += 1
return beta

X_b = np.hstack([np.ones((len(X), 1)), X])
initial_beta = np.zeros(X_b.shape[1])
self._beta = gradient_decent(X_b, y, initial_beta, eta, n_iters, eps)
self.coef_ = self._beta[1:]
self.interception_ = self._beta[0]
return self

def fit_sgd(self, X, y, n_iters=5, t0=5, t1=50):
assert X.shape[0] == y.shape[0]
assert n_iters >= 1

def dL_sgd(beta, X_b_i, y_i):
coef = np.exp(beta.dot(X_b_i)) / (1 + np.exp(beta.dot(X_b_i)))
return (-y_i + coef) * X_b_i

def gradient_decent_sgd(X_b, y, initial_beta, n_iters, t0, t1):

def learning_rate(t):
return t0 / (t1 + t)

beta = initial_beta
m = len(X_b)
for i in range(n_iters):
indexes = np.random.permutation(m)
X_b_new = X_b[indexes]
y_new = y[indexes]
for j in range(m):
gradient = dL_sgd(beta, X_b_new[j], y_new[j])
beta = beta - learning_rate(i * m + j) * gradient
return beta

X_b = np.hstack([np.ones((len(X), 1)), X])
initial_beta = np.zeros(X_b.shape[1])
self._beta = gradient_decent_sgd(X_b, y, initial_beta, n_iters, t0, t1)
self.coef_ = self._beta[1:]
self.interception_ = self._beta[0]
return self

def predict_probability(self, X_predict):
X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
return self._sigmoid(X_b.dot(self._beta))

def predict(self, X_predict):
assert self.coef_ is not None and self.interception_ is not None
assert X_predict.shape[1] == len(self.coef_)
probability = self.predict_probability(X_predict)
return np.array(probability >= 0.5, dtype='int')

def score(self, X_test, y_test):
y_predict = self.predict(X_test)
return np.sum(y_predict == y_test) / len(y_test)

def __repr__(self):
return "LogitRegression()"

测试代码如下

from LogitRegression import LogitRegression
from sklearn.model_selection import train_test_split
from sklearn import datasets
from sklearn import preprocessing

breast_cancer = datasets.load_breast_cancer()
y = breast_cancer.target
X = breast_cancer.data
train_X, test_X, train_y, test_y = train_test_split(X, y, test_size=0.3, random_state=123)

scaler = preprocessing.StandardScaler()
train_X = scaler.fit(train_X).transform(train_X)
test_X = scaler.transform(test_X)

logit = LogitRegression()
logit.fit(train_X, train_y)
print(logit.score(test_X, test_y))

运行结果如下所示：