统计学习方法 学习笔记--6.1.1 Logistic 回归 （梯度下降法）

统计学习方法 学习笔记--6.1.1  Logistic 回归 （梯度下降法）

目标：

已知训练数据 x1，x2， z。未来给你 x1 和 x2 ，判断 z 的值

西瓜密度	含糖率	好瓜与否	事件
x1	x2	z
0.697	0.46	1	z(1)
0.774	0.376	1	z(2)
0.634	0.264	1	z(3)
0.608	0.318	1	z(4)
0.556	0.215	1	z(5)
0.403	0.237	1	z(6)
0.481	0.149	1	z(7)
0.437	0.211	1	z(8)
0.666	0.091	0	z(9)
0.243	0.267	0	z(10)
0.245	0.057	0	z(11)
0.343	0.099	0	z(12)
0.639	0.161	0	z(13)
0.657	0.198	0	z(14)
0.36	0.37	0	z(15)
0.593	0.042	0	z(16)
0.719	0.103	0	z(17)

手段：

模型： z = w0+ x1*w1 + x2*w2

求出最佳的 w0, w1, w2

思考：

因为 z 的值是 0 和 1 二维的，可以考虑使用 sigmoid 函数： σ(z) = 1/(1+exp(-z))

当 z 的值为1的概率 如下：

P(Z =1 | x) = 1/ (1+exp-(w0+ x1*w1 + x2*w2) )

= exp(w0+ x1*w1 + x2*w2) / (1+ exp(w0+ x1*w1 + x2*w2))

P(Z=0|x) = 1 - P(Y=1|x) = 1/ (1+ exp(w0+ x1*w1 + x2*w2))

[说说自己的理解：

当 Z ->1时，σ(1) -> exp(z) / (1+exp(z));

当 Z ->0时，σ(0) -> exp(0) / (1+exp(z)) = 1 / exp(z); ]

模型参数估计：

Logistic 模型：

P(Z=1 | x) = π(x) = exp(w0+ x1*w1 + x2*w2) / (1+ exp(w0+ x1*w1 + x2*w2)) ,

P(Z=0 | x) = 1 - π(x) = 1/ (1+ exp(w0+ x1*w1 + x2*w2)) ,

似然函数： Π [ π(xi) ]^yi * [1- π(xi) ]^(1-yi)

[复习：似然函数定义：

函数的分布模型已经固定，当事件 z(1) ~ z(17)全部发生的概率 L(z1,z2,...,z17; w)，即联合概率密度。

L 越大，概率越大，模型的参数越合理 。

对似然函数取对数，新函数 L(w) 的最大值也为L的最大值]

对数似然函数为：

L(w) = ∑[yi*log( π(xi) )+(1-yi)*log(1- π(xi) )]

=∑[yi* (w ∙ xi) - log(1+exp (w ∙ xi) ]

问题转变：

问题由模型参数选择，变为选取合适的w使得L(w)最大化，即目标函数最优化的问题

运用梯度下降法：

对于线性函数 :

ΔL= Δw * L'(w)

由Taylor公式得到，对于其他函数：

ΔL <-- Δw * L'(w)

预定 初始值 w0，逐步迭代，每步移动 ∆w，当遇到极大值的时候 y′ (x) ->0, ∆y ->0

w(n) = w(n-1) + L' (w(n-1)) * ∆w

之前的理解 Δw是有问题的，新理解：

[

w n+1 = w n + λ * ∇f(w n ) ，其中 η 为步长，∇f(w n )为梯度

公式中的 wn ，事实上为元组，包含w1, w2...等多个维度

假设 wn =(w_x, w_y)，横轴为w_x，纵轴为w_y

垂直于纸面方向为 L(wx, wy)

则 wn+1 = wn + Δw

其中，Δw分为w_x 和 w_y两个分量，为保证 Δw 移动方向与等高线垂直。

则Δw_x = λ * ∂L/∂w_x, Δw_y = λ * ∂L/∂w_y

则 y 和 x 两个方向的 delta 值的比值为，Δw_y /Δw_x = ∂w_x/∂w_y，该方向与等高线的方向垂直，为最佳的迭代方向。

使用 偏导 ∂L/∂w_x 和 ∂L/∂w_y的目的，主要是协调 x 和 y两个方向上迭代的步调一致，实现总体方向最优。

再使用步长 λ，调整每次迭代跨过的距离。

在《统计学习方法》中：

一维搜索 λk，使得 L(w n + λ k *∇f(w n )) = min L(w n + λ*∇f(w n ))

从而获得最快的迭代速度，达到 L min

]

-----

故，当 w(n) 接近最优解时，

故，当 w(n) 接近最优解时，

L(w n) ≈ L(w n-1)

迭代终止条件为

L(w n) - L(w n-1) < 10^(-8)

计算 w(n) 迭代公式：

由 L(w) 求导

L(w)= ∑[zi (w ∙ xi )-log(1+exp(w ∙ xi ))]

导数：

∂L/∂w=∑[ xi *  [ zi - 1/(1+exp(-(w ∙x)) )] ]

带入 sigmoid函数：

∂L/∂w=∑[ xi *  [zi- σ(w∙x)] ]

-->

w(n) = w(n-1) + ∆w * X ∙ ( z - sigmoid(w ∙ x))

代码如下：

#西瓜书程序 watermelon

def loadDataSet_w():

dataMat = []; labelMat = []

fr = open('watermelon_1.txt') #原始数据testSet的数据格式为：数据1，数据2，label

for line in fr.readlines():

lineArr = line.strip().split()

dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) #其中的1.0为迭代的初始值

labelMat.append(int(lineArr[2]))

return dataMat,labelMat

dataArr_w, labelMat_w = loadDataSet_w()

weights_w = gradAscent(dataArr_w, labelMat_w)

def sigmoid(inX):

return 1.0/(1+exp(-inX))

#利用已有数据训练模型，目的是得到 weight向量 w = [w1, w2, ..., wn], 使得 dataMatIn * weights = classLabels

def gradAscent(dataMatIn, classLabels):

dataMatrix = mat(dataMatIn) #convert to numpy.matrixlib.defmatrix.matrix

labelMat = mat(classLabels).transpose() #convert to NumPy matrix，再转置

m,n = shape(dataMatrix)

alpha = 0.001 #alpha 步长

maxCycles = 500 #maxCycles 迭代次数

weights = ones((n,1)) #初始迭代，w向量赋值为 1

for k in range(maxCycles): #heavy on matrix operations，跑500次

h = sigmoid(dataMatrix*weights)

error = (labelMat - h)

weights = weights + alpha * dataMatrix.transpose()* error

return weights

dataArr_w, labelMat_w = loadDataSet_w()

weights_w = gradAscent(dataArr_w, labelMat_w)

def plotBestFit_w(weights):

import matplotlib.pyplot as plt

dataMat,labelMat=loadDataSet_w()

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(0.1, 0.9, 0.1)

y = (-weights[0]-weights[1]*x)/weights[2]

ax.plot(x, y)

plt.xlabel('X1'); plt.ylabel('X2');

plt.show()

plotBestFit_w(weights_w.getA())

---------------------

说明：

原始数据来源于 《机器学习》（西瓜书）

Python代码，参考《机器学习实战》