前言
之前对什么是逻辑回归,以及它的公式由来做了说明。在明确了该分类器的函数式后,那么最佳的回归系数是多少呢?这是值得思考的问题,本篇博客将会对这个问题进行探讨。
回顾逻辑回归公式
逻辑回归 = 线性回归 + sigmoid函数
线性回归: z = w*x + b
sigmoid函数:y = 1 1 + e − z {1}\over{1+e^{-z}} 1+e−z1
逻辑回归:y = 1 1 + e − w ∗ x + b {1}\over{1+e^{-w*x + b}} 1+e−w∗x+b1
对于sigmoid函数,其输入z = w 0 w_0 w0 x 0 x_0 x0 + w 1 w_1 w1 x 1 x_1 x1 + w 2 w_2 w2 x 2 x_2 x2 + … + w n w_n wn x n x_n xn,写成向量的形式为:z = W T W^T WTX,即将两个数值向量对应的元素相乘再求和得到z。这里的X就是样本数据,即分类器输入的数据;W就是需要求的最佳参数(系数)。分类器最终的分类效果好不好,精确率怎么样,都会受W的影响,所以要尽可能寻找最佳的W。寻找最佳W,会涉及到最优化理论的相关知识。
梯度上升法
这是一个最优化算法,这个算法基于的思想是:要找到某函数的最大值,最好的方法是沿着该函数的梯度方向探寻。如果梯度记作▽,那么函数f(x, y)的梯度有如下表示:
∇
f
(
x
,
y
)
=
(
∂
f
(
x
,
y
)
∂
x
,
∂
f
(
x
,
y
)
∂
y
)
\nabla f(x, y)=\left(\frac{\partial f(x, y)}{\partial x}, \frac{\partial f(x, y)}{\partial y}\right)
∇f(x,y)=(∂x∂f(x,y),∂y∂f(x,y)),这个梯度表示要沿着x的方向移动
∂
f
(
x
,
y
)
∂
x
{\partial f(x, y)}\over{\partial x}
∂x∂f(x,y),沿着y的方向移动
∂
f
(
x
,
y
)
∂
y
{\partial f(x, y)}\over{\partial y}
∂y∂f(x,y)。需要注意的是,函数f(x, y)在待计算点必须有定义且可微。
梯度上升算法迭代公式:
w
:
=
w
+
α
∇
w
f
(
w
)
w:=w+\alpha \nabla_{w} f(w)
w:=w+α∇wf(w)
与之相反的是梯度下降算法,公式由加号变成减号:
w
:
=
w
−
α
∇
w
f
(
w
)
w:=w-\alpha \nabla_{w} f(w)
w:=w−α∇wf(w),且用来求函数最小值。
python代码实现梯度上升找最佳参数
现在想要完成的任务是,在一些样本点上,利用梯度上升法找最佳回归系数,即找逻辑回归中最佳参数w。
- 加载数据集
数据集下载网盘链接: https://pan.baidu.com/s/1S9h2BMAyJejd6r28pQwprA
提取码: qp6y
数据集:这是一个包含100个样本点的数据集,每个点包含两个特征(数值型)。即是一个100X3(100行3列)的矩阵,第三列是样本标签(类别)。
def load_dataset():
"""加载数据集"""
dataSet = [] # 数据
labelSet = [] # 标签
fr = open('lr-testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
# 为了便于后续计算,将数据集第一列值设为1.0
dataSet.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelSet.append(int(lineArr[2]))
return dataSet, labelSet
- sigmoid函数
y = 1 1 + e − z {1}\over{1+e^{-z}} 1+e−z1
def sigmoid(X):
"""sigmoid函数"""
return 1.0/(1+exp(-X))
- 梯度上升法
这个函数有两个参数,第一个参数是100X3的矩阵,有100个样本,其中第一列 X 0 {X_0} X0是前面提到的便于后续计算,值设为1.0,后两列就是样本的2个特征 X 1 {X_1} X1和 X 2 {X_2} X2
def grad_ascent(dataSet, labelSet):
"""梯度上升算法"""
dataMatrix = mat(dataSet) # 转成numpy矩阵
labelMatrix = mat(labelSet).transpose() # 先转成numpy矩阵再进行矩阵转置
m, n = shape(dataMatrix)
print("样本数and特征数:", m, n)
alpha = 0.001 # 向目标移动的步长
maxCycles = 500 # 设置的最大的迭代次数
weights = ones((n, 1)) # 表示回归系数,3*1的矩阵,值全为1
for k in range(maxCycles):
h = sigmoid(dataMatrix*weights) # 矩阵相乘 # m*3 X 3*1
error = (labelMatrix - h) # m*1 - m*1
# alpha * dataMatrix.transpose() * error表示在每一列上的一个误差情况
# 得到x1,x2,xn的系数的偏移量
weights = weights + alpha * dataMatrix.transpose() * error
return weights
注意:ones函数可以用来创建任意维度和元素个数的数组,且元素值均为1;上面的代码中weights = ones((n, 1)) # 表示回归系数,3*1的矩阵,值全为1
表示我们自行创建了一个n*1的矩阵(n是特征数,为3,其中后两个是特征,第一个是设置的
X
0
X_0
X0)。因为我们最终实际所求的是回归系数w,因为有两个特征再加上设置的
X
0
X_0
X0,所以我们最后得到的是一个3X1的矩阵,即是对应的回归系数。
- 完整代码
from numpy import mat, shape, ones, exp # 可以用from numpy import * 代替此句
def load_dataset():
"""加载数据集"""
dataSet = [] # 数据
labelSet = [] # 标签
fr = open('lr-testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
if len(lineArr) == 1:
continue
# 为了便于后续计算,将数据集第一列值设为1.0
dataSet.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelSet.append(int(lineArr[2]))
return dataSet, labelSet
def sigmoid(inX):
"""sigmoid函数"""
return 1.0 / (1 + exp(-inX))
def grad_ascent(dataSet, labelSet):
"""梯度上升算法"""
dataMatrix = mat(dataSet) # 转成numpy矩阵
labelMatrix = mat(labelSet).transpose() # 先转成numpy矩阵再进行矩阵转置
m, n = shape(dataMatrix)
print("样本数and特征数:", m, n)
alpha = 0.001 # 向目标移动的步长
maxCycles = 500 # 设置的最大的迭代次数
weights = ones((n, 1)) # 表示回归系数,3*1的矩阵,值全为1
for k in range(maxCycles):
h = sigmoid(dataMatrix*weights) # 矩阵相乘 # m*3 X 3*1
error = (labelMatrix - h) # m*1 - m*1
# alpha * dataMatrix.transpose() * error表示在每一列上的一个误差情况
# 得到x1,x2,xn的系数的偏移量
weights = weights + alpha * dataMatrix.transpose() * error
return weights
# 测试
dataSet, labelSet = load_dataset()
lr_w = grad_ascent(dataSet, labelSet)
print("最佳回归系数:\n", lr_w)
输出结果:
样本数and特征数: 100 3
最佳回归系数:
[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
所以,现在求得的最佳的回归系数就是[ 0.48007329]和[-0.6168482 ]。
注意:
如果出现了如下错误信息:
return 1.0 / (1 + exp(-inX))
TypeError: only size-1 arrays can be converted to Python scalars
是因为在使用numpy数组的时候,我们最好也使用其中的数学函数处理数据。因为最开始在实现sigmoid函数的时候,exp我是从math模块中导入(from math import exp
)的,就出现了上边的错误信息。修改方法就是从numpy模块中导入这个exp。