假设一个数据集有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