Logistic 回归原理分析与 Python 实现

一、Logistic 回归原理分析

目标函数(sigmoid函数):

$$p(y=1|x;\theta)=h_\theta(x)=\frac{1}{1+\exp(\theta^Tx)}$$

这里的 $\theta$ 是模型参数,也就是回归系数,$\sigma$ 是sigmoid函数。这样 y=0 的“概率”就是:

$$p(y=0|x;\theta)=1-h_\theta(x)=\frac{exp(\theta^Tx)}{1+exp(\theta^Tx)}$$

对于逻辑斯蒂回归而言,可以得到如下的对数几率:

$$log(\frac{h_\theta(x)}{1-h_\theta(x)})=\theta_0+\theta_1x_1+\cdots+\theta_mx_m$$

在参数学习时,可以用极大似然估计方法求解。假设我们有n个独立的训练样本:

$$\{(x^{(1)}, y^{(1)}) ,(x^{(2)}, y^{(2)}), …, (x^{(n)}, y^{(n)})\}; y=\{0, 1\}$$

那每一个观察到的样本 $(x^{(i)}, y^{(i)})$ 出现的概率是:

$$p(y^{(i)},x^{(i)})=p(y^{(i)}=1|x^{(i)})^{y^{(i)}}(1-p(y^{(i)}=1|x^{(i)}))^{1-y^{(i)}}$$

对于整个样本集,每个样本的出现都是独立的,n个样本出现的似然函数为(n个样本的出现概率是他们各自的概率乘积):

$$L(\theta)=\prod p(y^{(i)},x^{(i)};\theta)=\prod p(y^{(i)}=1|x^{(i)};\theta)^{y^{(i)}}(1-p(y^{(i)}=1|x^{(i)};\theta))^{1-y^{(i)}}$$

对数似然函数为:

$$logit(\theta)=\sum (y^{(i)} logh_\theta(x^{(i)}) + (1-y^{(i)}) log(1-h_\theta(x^{(i)})) )$$

我们希望求得能使函数 $logit(θ)$ 取得最大值的参数 θ,因此我们设模型的代价函数(cost function)为:

$$l(\theta)=-logit(\theta)$$

这样,求解最佳参数 $\theta$,就可转化为求最小代价函数 $l(\theta)$。

二、步骤说明

1. 梯度下降法

1) 初始化参数 $\theta=(0,\theta_1,\theta_2,…,\theta_n )^T$,初始值为 $\theta=(0, 0, 0, …, 0 )^T$

2) 获取样本X,在每个样本前添加一个1,$x=(1,x_1,x_2,…,x_n )^T$

3) 设定步长α

4) 利用 sigmoid 函数对样本进行预测

5) 根据对数似然函数计算梯度 $g=\nabla_\theta l(\theta)$

6) 更新参数 $\theta=\theta-\alpha g$

重复步骤4~6,直到结果精度满足要求,或达到一定次数时迭代终止。

加入正则项的情况

代价函数增加正则项 $l(\theta)=-logit(\theta)+\frac{\lambda}{2m}\sum_{j=1}^m\theta_j^2$

因此步骤 5 修改为 $g=\nabla_\theta l(\theta)+\frac{\lambda}{m}\theta$

2. 牛顿法

1) 初始化参数 $\theta=(0,\theta_1,\theta_2,…,\theta_n )^T$,初始值为 $\theta=(0, 0, 0, …, 0 )^T$

2) 获取样本X,在每个样本前添加一个1,$x=(1,x_1,x_2,…,x_n )^T$

3) 利用 sigmoid 函数对样本进行预测

4) 根据对数似然函数计算梯度 $g=\nabla_\theta l(\theta)$

5) 根据对数似然函数计算 Hessian 矩阵H,$H_{ij}=\frac{\partial_2 l(\theta)}{\partial\theta_i \partial\theta_j}$

6) 更新参数 $\theta=\theta-H^{-1}g$

重复步骤 3~6,知道结果精度满足要求,或达到一定次数时迭代终止。

加入正则项的情况

修改步骤 4、5

4) $g=\nabla_\theta l(\theta)+\frac{\lambda}{m}\theta$

5) 计算 Hessian 矩阵 H,再让 H 加上下面这一项:

$$
\frac{\lambda}{m}
\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & ... & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
$$

三、代码实现

注:代码为 Python3 版本

完整的代码与样本放在了 Github 上,点此查看

1. 梯度下降

#############################################
# 梯度下降法估计参数
# X 是样本特征矩阵,每一行表示一个样本
# Y 是样本 X 中对应的标签,数组类型
# alpha 表示步长
#############################################
def gradientDescent(X, Y, alpha=0.001):
    # 在矩阵 X 的第一列插入 1
    X = np.insert(X, 0, 1, 1)
    # m 是训练样本数,n-1 是样本的特征数
    m, n = X.shape
    # 初始化 theta 值
    theta = np.mat(np.zeros(n))

    # 迭代求解theta
    for i in range(2000):
        theta -= alpha * gradient(X, Y, theta)

    debug("theta:", theta)
return theta

2. 梯度下降(带正则项)

#############################################
# 正则化的梯度下降法估计参数
# X 是样本特征矩阵,每一行表示一个样本
# Y 是样本 X 中对应的标签,数组类型
# alpha 表示步长
#############################################
def regularizedGradientDescent(X, Y, alpha=0.001):
    # 在矩阵 X 的第一列插入 1
    X = np.insert(X, 0, 1, 1)
    # m 是训练样本数,n-1 是样本的特征数
    m, n = X.shape
    # 初始化 theta 值
    theta = np.mat(np.zeros(n))
    # 惩罚系数
    lambda_ = 0.01 / Y.size

    # 迭代求解theta
    for i in range(200000):
        theta -= alpha * gradient(X, Y, theta)
        theta[0,1:] -= lambda_ * theta[0,1:]

    debug("theta:", theta)
return theta

3. 牛顿法

#############################################
# 牛顿法估计参数
# X 为特征矩阵,Y 为标签,iterNum 为迭代次数
#############################################
def newtonMethod(X, Y, iterNum=10):
    # 在矩阵 X 的第一列插入 1
    X = np.insert(X, 0, 1, 1)
    # m 是训练样本数,n-1 是样本的特征数
    m, n = X.shape
    # 初始化 theta 值
    theta = np.mat(np.zeros(n))

    # 迭代求解theta
    for iterIndex in range(iterNum):
        g = gradient(X, Y, theta)
        H = hessianMatrix(X, Y, theta)
        theta -= (H.I * g.T).T
        debug("theta({}):\n{}\n".format(iterIndex + 1, theta))

return theta

4. 牛顿法(带正则项)

#############################################
# 正则化牛顿法估计参数
# X 为特征矩阵,Y 为标签,iterNum 为迭代次数
#############################################
def regularizedNewtonMethod(X, Y, iterNum=10):
    # 在矩阵 X 的第一列插入 1
    X = np.insert(X, 0, 1, 1)
    # m 是训练样本数,n-1 是样本的特征数
    m, n = X.shape
    # 初始化 theta 值
    theta = np.mat(np.zeros(n))
    # 惩罚系数
    lambda_ = 0.01 / Y.size
    # hessian 矩阵正则项
    A = np.eye(theta.size)
    A[0,0] = 0

    # 迭代求解theta
    for iterIndex in range(iterNum):
        g = gradient(X, Y, theta) + lambda_ * theta
        H = hessianMatrix(X, Y, theta) + lambda_ * A
        theta -= (H.I * g.T).T
        debug("theta({}):\n{}\n".format(iterIndex + 1, theta))

return theta

四、测试结果

梯度下降法

牛顿法

相关文章

发表评论

电子邮件地址不会被公开。 必填项已用*标注