机器学习算法之： 逻辑回归 logistic regression (LR)

by joey周琦

LR介绍

逻辑回归属于probabilistic discriminative model这一类的分类算法

probabilistic discriminative mode这类算法的思路如下：

- 直接建模P(Ck|x)P(C_k|x)

- 利用最大似然估计和训练数据，估计出模型中的参数

该类想法相对于生成模型（probabilistic generated model) 有参数较少的优点。因为生成模型需要P(x|Ck)P(x|C_k)和先验概率P(Ck)P(C_k).

LR是工业界最长用的分类算法之一，其主要原因，个人认为有几点如下：

训练速度快，扛得住大数据

模型可解释度、可理解程度高，根据每个特征的系数，就可以判断出该特征在模型中的重要性，帮助判断模型是否合理

可以接受的精度

本文，对LR做一个简单的总结

LR二分类

首先做下简单的符号说明，在下述推导中，NN为样本个数，MM为特征数目，x∈RMx \in R_M为特征向量， KK为可分类的总数，P(Ck|x)P(C_k|x)表示在特征向量的前提下，判别为第k类的概率。w∈RMw \in R_M为LR模型的参数，即训练LR的主要目的就是要得到ww。

对于K=2K=2类，建模如下:

P(C1|x)=y(x)=σ(wTx)=11+e−wTxP(C_1|x)=y(x)=\sigma(w^Tx)=\frac{1}{1+e^{-w^Tx}}

σ\sigma是logistic sigmoid function, 它有个性质（1）：dσ(a)da=σ(1−σ)\frac{d\sigma(a)}{d a} = \sigma(1-\sigma)

其中yn=σ(an)y_n=\sigma(a_n)且 an=wTxa_n=w^Tx。假设我们有数据集{ xn,tnx_n,t_n }, n=1...Nn=1...N, NN为数据集的大小，nn表示第n个样本，tn=0or1t_n=0 or 1表示第nn个样本的真实类别，设t=(t1,…,tN)\mathbf{t}=(t_1,\dots,t_N),那么似然函数可以写为如下形式:

p(t|w)=∏n=1Nytnn(1−yn)1−tnp( \mathbf{t}| w)= \prod \limits_{n=1}^{N} y_n^{t_n}(1-y_n)^{1-t_n}

最大化似然函数等价于最大化对数似然函数可以写为

lnp(t|w)=∑n=1Ntnlnyn+(1−tn)ln(1−yn)\ln p( \mathbf{t}| w)= \sum \limits_{n=1}^{N} {t_n} \ln y_n + (1-t_n)\ln (1-y_n)

我们要最大化对数似然函数，就是等价于对数似然函数的相反数，则最小化目标函数可以写为：

J=−lnp(t|w)=−∑n=1Ntnlnyn+(1−tn)ln(1−yn)J = -\ln p( \mathbf{t}| w)= -\sum \limits_{n=1}^{N} {t_n} \ln y_n + (1-t_n)\ln (1-y_n)

目标函数对参数ww求导，利用上述σ\sigma函数的性质1，可以推到得到：

∂J∂w=∑n=1N(yn−tn)xn\frac{\partial J}{\partial w}=\sum \limits_{n=1}^{N} (y_n -t_n)x_n

可以利用梯度下降法对参数ww进行更新.更新公式为：

w=w−λ∂J∂w=w−λ∑n=1N(yn−tn)xnw=w-\lambda \frac{\partial J}{\partial w}=w-\lambda \sum \limits_{n=1}^{N} (y_n -t_n)x_n

其中λ\lambda是学习率，可以自己设置一个比较小的数字，或者通过其他算法计算（比如牛顿法）。可以看到每次更新都要有一个NN想得加和，这样很需要时间，所以更新策略可以采用“随机梯度下降法”，每次只用一个样本，如下：

w=w−λ∂J∂w=w−λ(yn−tn)xw=w-\lambda \frac{\partial J}{\partial w}=w-\lambda (y_n -t_n)x

这样大大节约了训练时间，也不会特别影响精度。有些系统为了避免过拟合，目标函数中可以加入L1或者 L2正则项（这里采用L2)，可以避免||w||2||w||^2过大,也就是为了避免模型为了拟合训练数据而变得过于复杂。加入正则项后，目标函数变为

J′=J+α2||w||2J' = J + \frac{\alpha}{2} ||w||^2

其中α\alpha为正则项惩罚系数。那么随机梯度的更新策略也就变为了

w=w−λ∂J′∂w=w−λ(yn−tn)xn−λαww=w-\lambda \frac{\partial J'}{\partial w}=w-\lambda (y_n -t_n)x_n-\lambda \alpha w

LR多分类

关于多分类的，在特征x∈RMx \in R_M下，模型判别为kk类的概率如下：

P(Ck|x)=yk(x)=exp(wTkx)∑Kj=1exp(wTjx)P(C_k|x)=y_k(x)=\frac{\exp(w_k^T x)}{ \sum \nolimits_{j=1}^K \exp(w_j^T x)}，

其中wj∈RMw_j \in R_M表示第j个模型的系数.另外，tn∈RK t_n \in R_K 表示的时第n个样本的真实分类, tnt_n向量只有一个元素为1，其余为0，若tnk=1t_{nk}=1则，表示该样本实际是第k类。T=t1,t2,…,tNT={t_1,t_2,\dots,t_N}表示所有样本的真实分类

那么最大似然函数则可以表示为：

P(T|w1,…,wK)=∏n=1N∏k=1KP(Ck|xn)tknP(T|w_1,\dots,w_K)= \prod_{n=1}^N \prod_{k=1}^K P(C_k | x_n)^{t_{kn}}

取负并取对数，则得到目标函数：

J=−logP(T|w1,…,wK)=−∑Nn=1∑Kk=1tnklogynkJ = -\log P(T|w_1,\dots,w_K) = - \sum \nolimits_{n=1}^N \sum \nolimits_{k=1}^K t_{nk} \log y_{nk}

可以推导（有点麻烦）：

∂J∂wj=∑Nn=1(ynj−tnj)xn \frac{\partial J }{\partial w_j} = \sum \nolimits_{n=1}^{N} (y_{nj} - t_{nj}) x_n

所以，梯度下降法对参数wjw_j进行更新公式：

wj=wj−λ∑Nn=1(ynj−tnj)xnw_j = w_j - \lambda \sum \nolimits_{n=1}^{N} (y_{nj} - t_{nj}) x_n

随机梯度下降法：wj=wj−λ(ynj−tnj)xnw_j = w_j - \lambda (y_{nj} - t_{nj}) x_n

随机梯度下降法+L2正则，目标函数变为：J′=J+α2||w||2J'=J + \frac{\alpha}{2} ||w||^2

则参数的更新变为：

wj=wj−λ(ynj−tnj)xn−αλwjw_j = w_j - \lambda (y_{nj} - t_{nj}) x_n - \alpha \lambda w_j

值得注意的是：

数据来自kaggle练习比赛的https://www.kaggle.com/c/digit-recognizer/data, 在这里我只用了train.csv数据，其中10%作为测试数据

由于数据的feature是像素，0~255，比较大，所以学习率要设置的低一些，或者将feature归一化。

加入正则项并不一定可以提高精度

结果如果2分类，区分数字0与9，正确率99%。如果是多分类，正确率为84%左右。

LR二分类

import numpy as np

p_test = [ ]
p_train = [ ]
lines = [line.rstrip('\n') for line in open('train.csv')]
label = 0
features = np.zeros(784) #28*28 feature

for line in lines[1:]:
line = line[0:-1]
data = line.split(",")
for i in range(len(data)):
data[i]=float(data[i])
data=np.array(data)
if np.random.rand()<0.1:
p_test.append(data)
else:
p_train.append(data)

#Train LR for 2 class
learing_r = 0.001
w = np.random.rand(784) - 0.5
for img in p_train:
label = img[0]
value = np.array(img[1:])
if label==9: # classify 9 and 0 , ignore others
t=1
elif label==0:
t=0
else:
continue
y = 1.0/(1+np.exp(-np.dot(w,value)))
w = w - learing_r * (y-t) * value

#Test LR
right = 0
wrong = 0
for img in p_test:
label = img[0]
value = np.array(img[1:])
if label==9: # classify 9 and 0 , ignore others
t=1
elif label==0:
t=0
else:
continue
y = 1.0/(1+np.exp(-np.dot(w,value)))
if y>0.5:
pt=1
else:
pt=0

if pt==t:
right+=1
else:
wrong+=1

print right*1.0/(right+wrong)

LR多分类

'''
Created on 2015 8 23

@author: joeyqzhou
'''
import numpy as np
import bigfloat

p_test = [ ]
p_train = [ ]
lines = [line.rstrip('\n') for line in open('train.csv')]
label = 0
features = np.zeros(784) #28*28 feature

for line in lines[1:]:
line = line[0:-1]
data = line.split(",")
for i in range(len(data)):
data[i]=float(data[i])
data=np.array(data)
if np.random.rand()<0.1:
p_test.append(data)
else:
p_train.append(data)

#Train LR for multiple class
learing_r = 0.000001
W = [ ]
K = 10 # 10 class
for i in range(K):
W.append(np.random.rand(784)*0.1 - 0.05)
for img in p_train:
label = img[0]
value = np.array(img[1:])
y = np.zeros(K)
temp_y = np.zeros(K) #unnormalized y
for i in range(K):
temp_y[i] = bigfloat.exp( np.dot(W[i], value) )
for i in range(K):
y[i] = temp_y[i]/np.sum(temp_y)
for i in range(K): #update the parameter
if abs(i-label) < 0.0001:
W[i] = W[i] - learing_r * (y[i]-1) * value
else:
W[i] = W[i] - learing_r * (y[i]-0) * value

#Test LR for multiple class
right_count = 0
whole_count = 0
for img in p_test:
label = img[0]
value = np.array(img[1:])
y = np.zeros(K)
for i in range(K):
y[i] =  np.dot(W[i], value)  #not necessary to normalize

inx = np.argmax(y)

if abs(inx - label)<0.001:
right_count += 1
whole_count += 1

print right_count*1.0/whole_count