上一篇:机器学习之python实现梯度下降
本文禁止转载!
网上许多教程都是使用python的一些包来使用线性回归,根本不能让我们搞清楚线性回归的原理。本文将使用纯python实现线性回归。
一、原理
线性回归是用一个一次函数来拟合。比如我们有100条奶茶销售量与当天气温的记录,我们画成散点图如下:
横坐标是当日气温,纵坐标是奶茶销售量。可以看到数据近似分布在一条直线上。于是我们假设它满足函数:
y
=
a
x
+
b
y=ax+b
y=ax+b
那么我们的目标就是求出一个适当的a,b。
我们随便蒙一个a=-0.6,b=16。画图如下:
好像还不错哈,但是这样主观估计肯定不行,我们要用数学方法找到最合适的解。
损失函数:
所以我们要有一个方法来评估一条直线与数据的拟合程度,这个方法我们称之为损失函数。
我们计算出每个点到直线的距离(y方向的距离),将它们加起来,我们将得到的数称为误差(loss),显然,误差越小,说明这条直线拟合程度越好。
我们设有n个数据点:(x1,y1),(x2,y2),…,(xi,yi),…,(xn,yn)。
将Yi记为第i个点的预测值,即:
Y
i
=
a
x
i
+
b
Y_i=ax_i+b
Yi=axi+b
所以:
l
o
s
s
=
∑
i
=
0
n
∣
Y
i
−
y
i
∣
loss=\sum_{i=0}^n{|Y_i-y_i|}
loss=i=0∑n∣Yi−yi∣
或
或
或
l
o
s
s
=
∑
i
=
0
n
(
Y
i
−
y
i
)
2
loss=\sum_{i=0}^n{(Y_i-y_i)^2}
loss=i=0∑n(Yi−yi)2
用后者比较好,毕竟绝对值不好展开。
所以:
l
o
s
s
=
∑
i
=
0
n
(
a
x
i
+
b
−
y
i
)
2
loss=\sum_{i=0}^n(ax_i+b-y_i)^2
loss=i=0∑n(axi+b−yi)2
l
o
s
s
=
∑
i
=
0
n
(
a
2
x
i
2
+
2
a
b
x
i
−
2
a
x
i
y
i
+
b
2
−
2
b
y
i
+
y
i
2
)
loss=\sum_{i=0}^n(a^2x_i^2+2abx_i-2ax_iy_i+b^2-2by_i+y_i^2)
loss=i=0∑n(a2xi2+2abxi−2axiyi+b2−2byi+yi2)
l
o
s
s
=
a
2
∑
i
=
0
n
x
i
2
+
2
a
b
∑
i
=
0
n
x
i
−
2
a
∑
i
=
0
n
x
i
y
i
+
n
b
2
−
2
b
∑
i
=
0
n
y
i
+
∑
i
=
0
n
y
i
2
loss=a^2\sum_{i=0}^n x_i^2 + 2ab\sum_{i=0}^n x_i -2a\sum_{i=0}^n x_iy_i + nb^2 -2b\sum_{i=0}^n y_i + \sum_{i=0}^ny_i^2
loss=a2i=0∑nxi2+2abi=0∑nxi−2ai=0∑nxiyi+nb2−2bi=0∑nyi+i=0∑nyi2
有没有感觉很熟悉,没错,他就是最小二乘法。
在这个式子中,a和b是未知的,所以它是关于a,b的一个函数。
而我们的目的就是求得一个a,b使得loss最小。
我们可以直接用数学方法(最小二乘法)计算出满足要求的a,b。
但是我们现在使用梯度下降的方法来找到a,b。
这就要用到我们上一篇文章所写的机器学习之python实现梯度下降。
二、代码实现
1.导入必要的包
这里我们不使用jupyter notebook,因为它画动图似乎有些问题,如果你想使用的话,参考:jupyter notebook matplotlib绘制动态图并显示在notebook中
import matplotlib.pyplot as plt
# %matplotlib inline
import numpy as np
import random
import shutil,os
2.生成数据并画出散点图
f = lambda x:-0.5*x+15
X = np.random.random(100)*50-15
Y = (np.array(list(map(f, X)))+np.random.normal(loc=0.0, scale=2.0, size=(100,))).astype('int32')
Y[Y<0]=0
plt.plot(X,Y,'.')
plt.show()
3.构建损失函数
xi_2 = np.sum(X**2)
yi_2 = np.sum(Y**2)
xiyi = np.sum(X*Y)
xi = np.sum(X)
yi = np.sum(Y)
n = len(X)
# 损失函数
loss = lambda a,b:xi_2*a**2 + n*b**2 + 2*xi*a*b - 2*xiyi*a - 2*yi*b + yi_2
# 损失函数在a,b方向上的导数
loss_da = lambda a,b:2*xi_2*a + 2*xi*b - 2*xiyi
loss_db = lambda a,b:2*n*b + 2*xi*a - 2*yi
4.梯度下降
# 设置学习率
learning_rate = 0.000025
# 精确度
precision = 0.001
# 随机生成一个初始a,b
min_a = 1
min_b = 0.01
# 计数
count = 0
# 记录过程
ab_list = []
while True:
ab_list.append([min_a,min_b])
# 求出该点导数
da_ = loss_da(min_a,min_b)
db_ = loss_db(min_a,min_b)
print('count:{},a:{},b:{}'.format(count,min_a,min_b))
# 沿导数方向前进
step_a = learning_rate * da_
step_b = learning_rate * db_ * 10
min_a -= step_a
min_b -= step_b
count += 1
if step_a**2+step_b**2 < precision**2 or count > 5000:
ab_list.append([min_a,min_b])
break
print('count:{},a:{},b:{}'.format(count,min_a,min_b))
5.画出动态图
img_path = './img/'
shutil.rmtree(img_path, ignore_errors=True)
os.makedirs(img_path)
# %matplotlib tk
ab_np = np.array(ab_list)
A = ab_np[:,0].reshape([-1])
B = ab_np[:,1].reshape([-1])
LINE_X = np.linspace(-15,35,100)
fig, ax = plt.subplots()
plt.xlim(-15, 35)
plt.ylim(-5, 25)
y1 = []
for i in range(len(A)):
ax.cla()
plt.xlim(-15, 35)
plt.ylim(-5, 25)
ax.plot(X,Y,'.')
ax.plot(LINE_X,A[i]*LINE_X+B[i])
plt.savefig("{}{:05d}.png".format(img_path,i))
plt.pause(0.02)
plt.show()
6.保存成gif
from PIL import Image
import imageio
filenames = os.listdir(img_path)
filenames.sort()
frames = []
for image_name in filenames:
image = Image.open(img_path+image_name)
array = np.array(image)
frames.append(array)
imageio.mimsave('output.gif', frames, 'GIF', duration=0.1)