假設一個數據集有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