一、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个独立的训练样本{(x1, y1) ,(x2, y2), …, (xn, yn)}; y={0, 1}。那每一个观察到的样本(xi, yi)出现的概率是：
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)})=\prod p(y^{(i)}=1|x^{(i)})^{y^{(i)}}(1-p(y^{(i)}=1|x^{(i)}))^{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) 初始化参数 \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
三、牛顿法
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
