python多线性回归模型 python线性回归模型图

网上许多教程都是使用python的一些包来使用线性回归，根本不能让我们搞清楚线性回归的原理。本文将使用纯python实现线性回归。

一、原理

线性回归是用一个一次函数来拟合。比如我们有100条奶茶销售量与当天气温的记录，我们画成散点图如下：

横坐标是当日气温，纵坐标是奶茶销售量。可以看到数据近似分布在一条直线上。于是我们假设它满足函数：

$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=ax_i+b$
所以：
$loss=\sum_{i=0}^n{|Y_i-y_i|}$
$或$
$loss=\sum_{i=0}^n{(Y_i-y_i)^2}$
用后者比较好，毕竟绝对值不好展开。
所以：
$loss=\sum_{i=0}^n(ax_i+b-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=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$
有没有感觉很熟悉，没错，他就是最小二乘法。
在这个式子中，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)