机器学习算法学习之线性回归
第一段是心路历程,不想看的可以跳到第二段。研究生开学以后自己的学习过程总是很急躁,缺乏细心沉淀,虽然又重新开始写博客了,不过对基础知识的学习总是有点儿操之过急,之前不想看机器学习算法,直接跳到了深度学习学习。最近由于实验室工作原因,又要学习机器学习方面的知识,学习一段时间后发现学习机器学习能积累很多理论基础的。那么现在换个方式,把基础打牢,慢慢沉淀吧。
废话说完了,本文会从线性回归谈起,会涉及很多数学内容,但是也会尽量通过一些函数图形让大家(和以后复习的自己)容易理解。随后才会正式进入逻辑回归,然后再尽可能来进行一些实战。当然了,自己也是萌新,希望有错误的地方大佬们提供一些指正。
初探线性回归
啥叫线性回归呢,我们从个通俗易懂的例子开始说起:
Living area( m 2 m^2 m2) | beadroom | Price(万RMB) |
---|---|---|
195 | 3 | 280 |
150 | 3 | 231 |
220 | 3 | 258 |
130 | 2 | 162 |
280 | 4 | 380 |
假设上面是一组住房面积、居室数和房屋总价的关系。房屋总价不光是由住房面积因素决定的,也与居室数、交通是否便利、基础设施是否完备等因素决定,我们这里只考虑两个因素(面积和居室数)。为了更可观查看数据,下图给出房屋面积和房屋总价在二维图像上的表示(为了方便画图图像相当于只考虑了面积因素):
如图,横坐标代表房屋面积,纵坐标代表房屋总价,这里只是大致在(130,162)和(280,380)两点之间画了一条直线,我们可以发现其他离散点也大致在这条直线的周围。
那么,现在根据上面的函数图像,在任意给出一个纵坐标(方面面积),就可以在这条直线上给出一个对应的纵坐标(房屋总价),这个纵坐标就是我们对面积为x的房屋总价的预测。二维的线性回归的方程大家在高中应该都学过:
y = a x + b y=ax+b y=ax+b
显然,如果我们考虑居室数的元素,二维的情况是不够的,例如我们把房屋面积、居室数分别用 x 1 x_1 x1, x 2 x_2 x2来表示,我们仍将其当做线性方程看待,那么我们可以把二维的方程推广到三维:
y = θ 0 + θ 1 x 1 + θ 2 x 2 y = \theta_0 + \theta_1x_1 + \theta_2x_2 y=θ0+θ1x1+θ2x2
现在,如果我们把它推广到多维:
h θ ( x ⃗ ) = [ θ 0 θ 1 ⋮ θ n ] T ⋅ [ 1 x 1 ⋮ x n ] h_{\theta}(\vec{x}) = \left[ \begin{matrix} \theta_0\\ \theta_1\\ \vdots\\ \theta_n\\ \end{matrix} \right]^T \cdot \left[ \begin{matrix} 1\\ x_1\\ \vdots\\ x_n\\ \end{matrix} \right] hθ(x)=⎣⎢⎢⎢⎡θ0θ1⋮θn⎦⎥⎥⎥⎤T⋅⎣⎢⎢⎢⎡1x1⋮xn⎦⎥⎥⎥⎤
现在我们把上面的公式简化来写,其中 θ ⃗ \vec{\theta} θ和 x ⃗ \vec{x} x都是列向量:
h θ ( x ⃗ ) = θ ⃗ T x ⃗ h_{\theta}(\vec{x}) = \vec{\theta}^T\vec{x} hθ(x)=θTx
Okay,现在我们已经有了基础的公式作为支撑,那么给定一堆离散点,如我们上面表格中的(195,3,280),(150,3,231)等离散点(可能会很多很多点),如何找到一个最优的线性方程,即找到合适的 θ ⃗ T \vec{\theta}^T θT,使得我们预测的效果更好呢?先清晰一下我们的目的,我们希望的所预测的点和我们的函数偏离度尽可能小,以我的理解,我们可以用概率学中的概念来理解,拿就是方差尽可能的小。在机器学习中,我们定义线性回归的损失函数(loss function)如下:
J ( θ ⃗ ) = 1 2 ∑ i = 1 m ( h θ ⃗ ( x ( i ) ) − y ( i ) ) 2 J(\vec{\theta}) = {1\over2}{\sum^{m}_{i=1}(h_{\vec{\theta}}(x^{(i)})-y^{(i)})^2} J(θ)=21i=1∑m(hθ(x(i))−y(i))2
其中 x ( i ) x^{(i)} x(i)是样本的第i个样本的属性(注意 x x x是是个列向量,其中可能包括多种属性,如我们之前例子中可以包括房屋面积、居室数等属性), y ( i ) y^{(i)} y(i)是指第个样本对应的值(如上例中指房价)。这个公式可能猛地一看很懵逼,实际上我们仔细观察,它和样本的方差很类似,为了方便大家回顾,样本的方差公式为:
S = 1 n − 1 ∑ i = 1 n ( X ‾ − X i ) 2 S = {1\over{n-1}}{\sum^{n}_{i=1}(\overline{X}-X_i)^2} S=n−11i=1∑n(X−Xi)2
之所以前面的分数换成了 1 2 1\over2 21,其实可以在求导之后更容易消掉。
我们知道 y = x 2 y=x^2 y=x2类型的函数是凸函数,其函数图像大致为以下的样子:
凸函数有个重要的属性就是二阶导横大于0,也就是说,当一阶导为0时,就是其极小值。二阶指数函数的二阶导在定义域上横大于零,因此此时求到的极小值就是其全局最小值。通过上面 y = x 2 y=x^2 y=x2的函数图像我们也可以清晰的理解这个道理。
回到之前的内容,现在我们的目标是尽可能使我们预测的准确。就是在样本给定的情况下求 J ( θ ) J(\theta) J(θ)的最小值,将上面一阶导为0的情况推广到多维,就是使关于 θ \theta θ的函数 J ( θ ) J(\theta) J(θ)梯度为0。(注意 θ \theta θ是个列向量,即包含 θ 0 \theta_0 θ0、 θ 1 \theta_1 θ1…… θ n \theta_n θn,梯度为0时 ∂ h ∂ θ 0 = 0 {\partial{h}\over\partial{\theta_0}}=0 ∂θ0∂h=0、 ∂ h ∂ θ 1 = 0 {\partial{h}\over\partial{\theta_1}}=0 ∂θ1∂h=0……、 ∂ h ∂ θ n = 0 {\partial{h}\over\partial{\theta_n}}=0 ∂θn∂h=0)
即:
X ⃗ T X ⃗ θ ⃗ = X ⃗ y ⃗ \vec{X}^T\vec{X}\vec{\theta} = \vec{X}\vec{y} XTXθ=Xy
如果上面的 X ⃗ T X ⃗ \vec{X}^T\vec{X} XTX可逆,那么可以算出:
θ ⃗ = ( X ⃗ T X ⃗ ) − 1 X ⃗ y ⃗ \vec{\theta}=(\vec{X}^T\vec{X})^{-1}\vec{X}\vec{y} θ=(XTX)−1Xy
如果不可逆,我们就可以加一个很小的正数偏置 λ \lambda λ,使其可逆:
θ ⃗ = ( X ⃗ T X ⃗ + λ E ) − 1 X ⃗ y ⃗ \vec{\theta}=(\vec{X}^T\vec{X}+\lambda{E})^{-1}\vec{X}\vec{y} θ=(XTX+λE)−1Xy
这里稍微解释一下为啥加了 λ \lambda λ就可逆了,因为其实 X ⃗ T X ⃗ \vec{X}^T\vec{X} XTX是半正定的,如我们在两边乘上 μ ⃗ \vec{\mu} μ:
μ ⃗ T X ⃗ T X ⃗ μ ⃗ = ( X ⃗ μ ⃗ ) T X ⃗ μ ⃗ = μ 1 2 + μ 2 2 + … … + μ n 2 ≥ 0 \vec{\mu}^T\vec{X}^T\vec{X}\vec{\mu} = (\vec{X}\vec{\mu})^T\vec{X}\vec{\mu}=\mu_1^2+\mu_2^2+……+\mu_n^2\geq0 μTXTXμ=(Xμ)TXμ=μ12+μ22+……+μn2≥0
如果不可逆,对于一个不等于0向量的 μ ⃗ \vec\mu μ可以使上式存在等于0的情况,我们只需加一个小偏置 λ \lambda λ那么不等于0的情况就不存在了。
回顾
从上面的一个很简单的例子开始,我们一步步理解了线性回归的含义,并大致理解了如何寻找合适的
θ
⃗
\vec{\theta}
θ来使我们的预测模型尽可能的准确,但是其实我们面临一个难题:如何尽可能高效的求矩阵的逆?
对于阶数较低的情况,用程序求逆很简单,但是如果样本很大,我们上面
X
⃗
\vec{X}
X的阶数就会随之增加。因此我们需要一个更快、更高效的算法来计算
J
(
θ
⃗
)
J(\vec{\theta})
J(θ)的最小值,也就是我们接下来将介绍的梯度下降算法。
梯度下降
为了更好的理解算法,大家可以看下面的图:
上面仍是 y = x 2 y=x^2 y=x2的图像(这个函数比较好理解)。我们以上述图像为例来描述梯度下降的过程:假设我们有个绿色的小球一开始位于(-3,9)的位置,为了使其到达最低点,我们先使顺着图像滑动一部分,到x=-2.4左右的位置再继续下降,一直大约到x=0为止。这个过程就是所谓梯度下降的过程。
如何进行滑动呢?根据微积分我们知道,对于一个函数C,我们向v1方向移动 Δ v 1 \Delta{v1} Δv1,向v2方向移动 Δ v 2 \Delta{v2} Δv2,若两个变化量都很小时(注意重点!!!如果很大的话就不符合了,也就可能会出现梯度无限下降到负无穷的情况,当然也不是变化量约小越好,具体情况具体分析),那么:
Δ C ≈ ∂ C ∂ v 1 Δ v 1 + ∂ C ∂ v 2 Δ v 2 \Delta{C} \approx {{\partial{C}}\over{\partial{v1}}}\Delta{v1} + {{\partial{C}}\over{\partial{v2}}}\Delta{v2} ΔC≈∂v1∂CΔv1+∂v2∂CΔv2
我们又知道,C的梯度 ∇ C \nabla{C} ∇C为 ( ∂ C ∂ v 1 , ∂ C ∂ v 2 ) ({{\partial{C}}\over{\partial{v1}}},{{\partial{C}}\over{\partial{v2}}}) (∂v1∂C,∂v2∂C),因此我们可以把上式写成:
Δ C = ∇ C ⋅ v ⃗ \Delta{C} = \nabla{C}\cdot\vec{v} ΔC=∇C⋅v
如果我们取 v ⃗ = − η ∇ C \vec{v}=-\eta\nabla{C} v=−η∇C,那么:
Δ C = − η ∇ C ⋅ ∇ C = − η ∣ ∣ ∇ C ∣ ∣ \Delta{C} = -\eta\nabla{C}\cdot\nabla{C}=-\eta||\nabla{C}|| ΔC=−η∇C⋅∇C=−η∣∣∇C∣∣
其中我们知道 ∣ ∣ ∇ C ∣ ∣ ||\nabla{C}|| ∣∣∇C∣∣必然大于等于0,因此,我们找到了一个变化量 v ⃗ \vec{v} v,使我们滑动过程必然是下降的。这也就是梯度下降的原理。现在的问题就是如何求 ∇ C \nabla{C} ∇C了,我们可以对第i个 θ \theta θ求偏导数:
∂ ∂ θ i J ( θ ) = ( h θ ( x ⃗ ) − y ) x i ( 推 导 过 程 省 略 , 和 前 面 的 类 似 ) {{\partial}\over{\partial{\theta_i}}}J(\theta) = (h_{\theta}(\vec{x})-y)x_i (推导过程省略,和前面的类似) ∂θi∂J(θ)=(hθ(x)−y)xi(推导过程省略,和前面的类似)
其中 h θ ( x ⃗ ) h_{\theta}(\vec{x}) hθ(x)就是我们当前预测的结果,y是正确的预测结果。
算法实现
我从某二手房网站爬取了1000余条数据用于实战检测,当然,因为现实情况是,居室数和面积对房价只是部分影响,所以不一定完全符合线性模型,相关性大概训练完是87%,只是纯属娱乐!另外也测试了自己按线性模型生成的序列和机器学习实战中给出的符合线性模型的一组数据,这个效果当然比较好些,相关性可以达到99%,我们直接贴出全部代码,再逐步分析各个模块吧。我把代码放到github和相关训练集上了:
https://github.com/UnnameBao/machine_learing_code
#enconding:utf-8
'''
@Author: b0ring
@MySite: https://unnamebao.github.io/
@Date: 2019-11-13 15:47:04
@Version: 1.0.0
'''
import numpy as np
import random
#实现线性回归的类
class linear():
#在init中随机生成三个theta,并设置学习速率
def __init__(self):
self.theta = np.random.randn(3, 1)
#注意这个学习速率要根据情况变!因为房价的数值比较大,所以这里的theta很小
self.eta = 0.0000001
#epochs是训练次数,默认设置为1
def train(self,train_data,mini_batch_size,epochs=1):
for i in range(epochs):
#每次训练打乱train_data
random.shuffle(train_data)
#分批次训练,每次训练mini_batch_size个数据
mini_batches = [
train_data[k:k+mini_batch_size]
for k in range(0, len(train_data), mini_batch_size)]
for mini_batche in mini_batches:
#每次update进行梯度下降
self.update_mini_batch(np.array(mini_batche),self.eta)
#预测值,将theta和x做点乘
def predict(self,x):
return np.dot(x,self.theta)
#梯度下降算法实现
def update_mini_batch(self,mini_batch,eta):
#前两列是样本x的值
y = mini_batch[:,2:]
#theta0是常数,因此第一列需是1
new_col = np.array([[1]]*mini_batch.shape[0])
x = np.c_[new_col,mini_batch[:,:2]]
#计算当前预测值
y_hat = self.predict(x)
#对应公式求导公式,并更新新的theta
part_one = y - y_hat
self.theta = self.theta + (1/len(mini_batch))*np.sum((eta*part_one*x).T,axis=1).reshape((-1,1))
linear_ = linear()
f = open("data.txt","r")
lines = f.readlines()
data = [[float(one) for one in line.strip().split("\t")] for line in lines]
train_data = data[:15000]
test_data = data[15000:]
linear_.train(train_data,100)
input_ = np.array(data)[:,:2]
y = np.array(data)[:,2:]
new_col = np.array([[1]]*input_.shape[0])
input_ = np.c_[new_col,input_]
y_hat = linear_.predict(input_)
#计算预测值与真实值的相关系数
print(np.corrcoef(y_hat.T, y.T))
其他代码查看注释还是很清楚的,我们主要来看梯度下降的代码:
#梯度下降算法实现
def update_mini_batch(self,mini_batch,eta):
#前两列是样本x的值
y = mini_batch[:,2:]
#theta0是常数,因此第一列需是1
new_col = np.array([[1]]*mini_batch.shape[0])
x = np.c_[new_col,mini_batch[:,:2]]
#计算当前预测值
y_hat = self.predict(x)
#对应公式求导公式,并更新新的theta
part_one = y - y_hat
self.theta = self.theta + (1/len(mini_batch))*np.sum((eta*part_one*x).T,axis=1).reshape((-1,1))
上述代码对应我们之前推导出的公式:
∂
∂
θ
i
J
(
θ
)
=
(
h
θ
(
x
⃗
)
−
y
)
x
i
{{\partial}\over{\partial{\theta_i}}}J(\theta) = (h_{\theta}(\vec{x})-y)x_i
∂θi∂J(θ)=(hθ(x)−y)xi
Δ
C
=
−
η
∇
C
⋅
∇
C
\Delta{C} = -\eta\nabla{C}\cdot\nabla{C}
ΔC=−η∇C⋅∇C
由于theta0是常数项,因此我们在点乘之前需要在第一列设置为常数1。y_hat是我们当前模型对样本的预测值,y是我们样本的真实值,二者相减就对应公式中的 h θ ( x ⃗ ) − y h_{\theta}(\vec{x})-y hθ(x)−y。代码中因为我把学习速率 η \eta η设置的是正值,因此代码中我写的是y-y_hat,其实是一样的。
需要注意的是,因为我们通过分批次训练来修改theta值,所有有每次计算完 ( h θ ( x ⃗ ) − y ) x i (h_{\theta}(\vec{x})-y)x_i (hθ(x)−y)xi后要相加求均值(即我代码中每次训50个后计算平均下降的值)。
对于房价预测的相关系数:
如上所示,大约是87.4%(注意可能每次运行不一样,因为我theta初始是随机的,不过肯定不会偏差很多)。对应修改上面一部分代码,可以检验我通过线性模型随机生成的数据:
f = open("house.txt","r")
lines = f.readlines()
data = [[float(one) for one in line.strip().split("\t")] for line in lines]
train_data = data[:1000]
test_data = data[1000:]
预测线性模型的数据是很可观的:
后记
看原理就花时间蛮久的,总结时间好像更久,特别代码实现的时候遇到了一些小问题,主要一方面是使用numpy不太熟,另一方面是实践少,最开始预测房价的时候我碰到了无限下降的问题,因为自己把 η \eta η设置的值为0.01,对房价来说太大了,因为房价动辄几十万起步,明显就不是个很小的值。所以还要多动手,多实践,多交流!感谢你看到这里!