假设一个数据集有n个样本，每个样本有m个特征，样本标签y为{0, 1}。
数据集可表示为：
其中，x(ij)为第i个样本的第j个特征值，y(i)为第i个样本的标签。
X矩阵左侧的1相当于回归方程的常数项。
每个特征有一个权重（或系数），权重矩阵为：
开始可以将权重均初始化为1。
将特征及权重分别相乘得到Xw (即特征的线性组合，为n维列向量)。
经过Sigmoid函数处理得到预测值：
y为预测值（取值范围0-1），为n维列向量。
对于一个样本i，y(i)取值为1和0的概率分别是：
其中x(i)为第i个样本的特征向量，为m+1维行向量。
为了学习得到最佳的权重矩阵w，需要定义损失函数来优化。一个直观的想法是使用预测值y与观测值Y之间的误差平方和，但是这个损失函数是非凸函数，用梯度下降法不能得到全局极小值。所以我们采用最大似然法。
对于每一个样本，出现的概率为：
假设n个样本相互独立，概率相乘。似然函数为：
取对数，变乘法为加法，得到对数似然函数：
这就是我们需要最大化的目标函数。
梯度法
如采用梯度法，首先要对w求导：
其中，σ为Sigmoid函数。
最后使用梯度上升来更新权重：
其中α为步长。经过多次迭代后，求得似然函数的最大值及相应的w。
牛顿法
如采用牛顿法，需要计算二阶导数：
这是一个m×m的矩阵，称为Hessian矩阵，用H表示。
如果定义：
则：
根据牛顿迭代公式：
经过有限次迭代，达到收敛。
预测分类
如果用来预测分类，进行如下运算：
如y(i) > 0.5 判定为1，如y(i) < 0.5，判定为0
权重系数与OR的关系
下面讨论一下权重w与OR的关系。
根据OR的定义:
当其他特征值不变的情况下，某x(i)增加1，相应的和xw增加w(i)，OR值变为原来的exp(w(i)) 倍。
Python程序代码
from numpy import *
import matplotlib.pyplot as plt
# 加载数据def loadDataSet():
dataMat = []
labelMat = []
fr = open('data/testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split(',')
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
fr.close()
return dataMat, labelMat
# Sigmoid函数，注意是矩阵运算def sigmoid(inX):
return 1.0/(1+exp(-inX))
# 梯度上升算法def gradAscent(dataMatIn, classLabels):
dataMat = mat(dataMatIn)
labelMat = mat(classLabels).transpose()
m,n=shape(dataMat)
alpha = 0.01
maxCycles = 500
weights = mat(ones((n,1)))
weightsHis = [mat(ones((n,1)))] # 权重的记录，主要用于画图
for k in range(maxCycles):
h = sigmoid(dataMat*weights)
error = labelMat - h
weights = weights + alpha*dataMat.transpose()*error
weightsHis.append(weights)
return weights,weightsHis
# 简单的随机梯度上升，即一次处理一个样本def stocGradAscent0(dataMatIn, classLabels):
dataMat = mat(dataMatIn)
labelMat = mat(classLabels).transpose()
m,n=shape(dataMat)
alpha = 0.01
weights = mat(ones((n,1)))
weightsHis = [mat(ones((n,1)))] # 权重的记录，主要用于画图
for i in range(m):
h = sigmoid(dataMat[i]*weights)
error = labelMat[i] - h
weights = weights + alpha* dataMat[i].transpose() * error
weightsHis.append(weights)
return weights,weightsHis
# 改进的随机梯度算法def stocGradAscent1(dataMatIn, classLabels, numIter):
dataMat = mat(dataMatIn)
labelMat = mat(classLabels).transpose()
m,n=shape(dataMat)
alpha = 0.001
weights = mat(ones((n,1)))
weightsHis = [mat(ones((n,1)))] # 权重的记录，主要用于画图
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.001 # 动态调整alpha
randIndex = int(random.uniform(0,len(dataIndex))) # 随机选择样本
h = sigmoid(dataMat[randIndex]*weights)
error = labelMat[randIndex]- h
weights=weights + alpha * dataMat[randIndex].transpose() * error
del(dataIndex[randIndex])
weightsHis.append(weights)
return weights,weightsHis
# 牛顿法def newton(dataMatIn, classLabels, numIter):
dataMat = mat(dataMatIn)
labelMat = mat(classLabels).transpose()
m,n=shape(dataMat)
# 对于牛顿法，如果权重初始值设定为1，会出现Hessian矩阵奇异的情况.
# 原因未知，谁能告诉我
# 所以这里初始化为0.01
weights = mat(ones((n,1)))-0.99 weightsHis = [mat(ones((n,1))-0.99)] # 权重的记录，主要用于画图
for _ in range(numIter):
A = eye(m)
for i in range(m):
h = sigmoid(dataMat[i]*weights)
hh = h[0,0]
A[i,i] = hh*(1-hh)
error = labelMat - sigmoid(dataMat*weights)
H = dataMat.transpose() * A * dataMat # Hessian矩阵
weights = weights + H**-1 * dataMat.transpose() * error
weightsHis.append(weights)
return weights,weightsHis
def plotWeights(w): w = array(w) def f1(x):
return w[x,0,0]
def f2(x):
return w[x,1,0]
def f3(x):
return w[x,2,0] k = len(w)
x = range(0,k,1)
plt.plot(x,f1(x),'',x,f2(x),'',x,f3(x),'')
plt.show()
# 画出分类边界def plotBestFit(wei):
weights = wei.getA()
dataMat, labelMat = loadDataSet()
dataArr = array(dataMat)
n = shape(dataArr)[0]
xcord1=[]
ycord1=[]
xcord2=[]
ycord2=[]
for i in range(n):
if int(labelMat[i])==1:
xcord1.append(dataArr[i,1])
ycord1.append(dataArr[i,2])
else:
xcord2.append(dataArr[i,1])
ycord2.append(dataArr[i,2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = arange(-3.0,3.0,0.1)
y=(-weights[0]-weights[1]*x)/weights[2]
ax.plot(x,y)
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()
# 测试 data, labels = loadDataSet()
#weights,weightsHis = gradAscent(data, labels)#weights0, weightsHis0 = stocGradAscent0(data, labels)#weights1, weightsHis1 = stocGradAscent1(data, labels, 500)weights3, weightsHis3 = newton(data, labels, 10)
plotBestFit(weights3)print(weights3)
plotWeights(weightsHis3) |
运行结果：
1、梯度法迭代500次的分类边界及权重收敛情况
2、随机梯度法迭代500次的分类边界及权重收敛情况
2、随机梯度法迭代500次的分类边界及权重收敛情况
3、牛顿法迭代10次的分类边界及权重收敛情况，可以牛顿法要快很多。
转载于：http://blog.sina.com.cn/s/blog_44befaf60102wbbr.html