本文主要从数学原理入手,探讨如何通过已知数据,使电脑学习出线性回归方程。在上来开始详细介绍之前,我们先来看一个例子。
这是一个关于银行贷款额度的数据统计表。可以看出,在本例中银行贷款额度与贷款者的工资、年龄有关,换做计算机语言来讲,就是我们的预测函数受到两个变量的控制。如何才能得出本例的线性回归方程呢,首先从数学方面出发,通过从预测函数的函数式着手,对变量进行迭代,以下是详细步骤。
目录
一、损失函数的数学推导:
以下为损失函数的数学推导,不是必要内容,可以选择性跳过
在开始推导前,我们可以大致估算出我们的函数模型,函数模型的估算因人而异。不同的函数模型会影响最后的拟合效果,倘若函数模型过于复杂或者维度过高,则会出现过拟合的现象,本文不再予以讨论,如果读者感兴趣可以浏览作者有关正则化的博文,里面有详细的讲解。
实际上,引入矩阵的概念可以极大程度上方便我们的计算,如果不会矩阵的操作可以跳过这一点。通过矩阵的概念显然我们可以得出:预测函数=theta的转置 乘以 矩阵x
这里运用到了概率论的知识。我们得到了最后一行的结果,目的是要让P最大,显然需要让后面的式子越小越好,因而得出损失函数的数学概念:
对于我们是否带入矩阵操作,得到两个不同的表达式,这两个表达式都可以完成对线性回归方程的计算。
(1) 梯度下降法求线性回归方程
需注意,在梯度下降算法中,我们需要定义每次迭代的步长alpha,alpha的大小会影响迭代次数和是否能够找到损失函数的最值,通常取小于1的数。梯度下降算法具有一定的泛用性,在参数极多的情况下也能正常运行。
(2)正规方程法求线性回归方程
在正规方程法当中,我们无需去定义步长alpha的大小,也无需迭代,即可求得参数theta。但相反,正规方程需要求解参数x矩阵的转置,即便numpy函数库可以很方便的计算转置矩阵,但在参数极多的情况下,转置矩阵的运算量会达到一个十分恐怖的数量级,因而不具有泛用性。
二、源码展示
在梯度下降算法中,参数的数量仅涉及到对参数求导、参数迭代的过程,原理上类似。作者是通过notebook环境下编写的代码,下面展示源码。
导入必须的python库,我们用到了pandas、numpy和matplotlib库。
import pandas as pd #pandas库用来从csv表格中读取数据
import numpy as np
import matplotlib.pyplot as plt
定义求损失函数的函数模块,返回损失函数的函数值。
#求损失函数
def Loss(theta0, theta1, theta2, x1, x2, y):
h = theta0 + theta1 * x1 + theta2 * x2
diff = h - y #此时h和y都是一个5维矩阵,可以直接做减法,得到每一项与
原始值的差值
L = (diff.sum())**2 / (2 * x1.shape[0]) #x1.shape[0]得到的是矩阵的行数,在5维矩阵中即
表示有几个元素
return L
定义求各个参数偏导的函数模块,为梯度下降做准备。
#求损失函数中theta0的偏导
def partial_loss_theta0(theta0, theta1, theta2, x1, x2, y):
h = theta0 + theta1 * x1 + theta2 * x2
diff = h - y #求导的过程用上文推导的式子即可
partial = diff.sum() / x1.shape[0]
return partial
#求损失函数中theta1的偏导
def partial_loss_theta1(theta0, theta1, theta2, x1, x2, y):
h = theta0 + theta1 * x1 + theta2 * x2
diff = (h - y) * x1 #求导的过程用上文推导的式子即可
partial = diff.sum() / x1.shape[0]
return partial
#求损失函数中theta2的偏导
def partial_loss_theta2(theta0, theta1, theta2, x1, x2, y):
h = theta0 + theta1 * x1 + theta2 * x2
diff = (h - y) * x2 #求导的过程用上文推导的式子即可
partial = diff.sum() / x1.shape[0]
return partial
梯度下降核心算法如下。
#梯度下降求损失函数最小值点
def gradient_descent(x1, x2, y, alpha=0.1, theta0 = 1, theta1 = 1, theta2 = 1, steps = 5000):
num = 0 #迭代次数
result=[]
for step in range(steps):
num += 1
update0 = partial_loss_theta0(theta0, theta1, theta2, x1, x2, y)
update1 = partial_loss_theta1(theta0, theta1, theta2, x1, x2, y)
update2 = partial_loss_theta2(theta0, theta1, theta2, x1, x2, y)
theta0 = theta0 - alpha * update0 #对三个参数同时进行迭代
theta1 = theta1 - alpha * update1
theta2 = theta2 - alpha * update2
e = Loss(theta0, theta1, theta2, x1, x2, y)
result.append(e)
if e < 0.001:
break
return theta0,theta1,theta2,result
此处作者用到了对数据的特征缩放,目的是让数据更集中、更规则,从而大大提高了运算速度。此处附上notebook上的数据标准化后的结果展示:
ranges = pd.read_csv("ranges.csv")
ranges.Salary = (ranges.Salary - ranges.Salary.mean()) / ranges.Salary.std() #此处对原本的数据进行了标准化
ranges.Age = (ranges.Age - ranges.Age.mean()) / ranges.Age.std() #标准化后的数据会大大减少运算次数
ranges.Range = (ranges.Range - ranges.Range.mean()) / ranges.Range.std()
ranges #数据标准化后的结果展示
最后定义main函数,并通过matplotlib画出损失函数梯度下降的图像.
if __name__ == "__main__":
theta0,theta1,theta2,result = gradient_descent(ranges.Salary, ranges.Age, ranges.Range)
print("\n线性规划方程为:\n")
print("h(x) = %.3f + %.3f * x1 + %.3f * x2\n"%(theta0,theta1,theta2))
n = len(result)
x = range(n)
plt.plot(x,result,color='r',linewidth=3)
plt.title("Loss")
plt.xlabel("num")
plt.ylabel("y")
plt.show()
我们可以看到,经过标准化后的数据,经过大概三十次左右的的计算即可基本拟合出较好的预测函数,效率很高。