05 模型之母：简单线性回归的代码实现

2019-12-15 本文已影响0人 Japson

本文为饼干Japson原创，首发于公众号《数据科学家联盟》。更多大数据、机器学习、深度学习相关内容，敬请关注公众号。

0 前言

在《模型之母：简单线性回归&最小二乘法》中，我们从数学的角度理解了简单线性回归，并且推导了最小二乘法。

本文内容完全承接于上一篇，我们来以代码的方式，实现简单线性回归。话不多说，码起来

1 简单线性回归算法的实现

首先我们自己构造一组数据，然后画图

import numpy as np
import matplotlib.pyplot as plt

x = np.array([1.,2.,3.,4.,5.])
y = np.array([1.,3.,2.,3.,5,])

plt.scatter(x,y)
plt.axis([0,6,0,6])
plt.show()

15724267041605.jpg

下面我们就可以根据样本真实值，来进行预测。

实际上，我们是假设线性关系为： $f(x) = ax + b$ 这根直线，然后再根据最小二乘法算a、b的值。我们还可以假设为二次函数： $f(x) = ax^2 + bx + c$ 。可以通过最小二乘法算出a、b、c

实际上，同一组数据，选择不同的f(x)，即模型，通过最小二乘法可以得到不一样的拟合曲线。

不同的数据，更可以选择不同的函数，通过最小二乘法可以得到不一样的拟合曲线。

下面让我们回到简单线性回归。我们直接假设是一条直线，模型是： $\sum_{i=1}^{m} (y^{(i)}-\hat{y}^{(i)})^{2}$

根据最小二乘法推导求出a、b的表达式：

$a = \frac {\sum_{i=1}^{m} (x^{i}-\bar{x})} {\sum_{i=1}^{m} (x^{(i)}-\bar{x})^{2}} \qquad b=\bar{y} - a\bar{x}$

下面我们用代码计算a、b：

# 首先要计算x和y的均值
x_mean = np.mean(x)
y_mean = np.mean(y)

# a的分子num、分母d
num = 0.0
d = 0.0
for x_i,y_i in zip(x,y):   # zip函数打包成[(x_i,y_i)...]的形式
    num = num + (x_i - x_mean) * (y_i - y_mean)
    d = d + (x_i - x_mean) ** 2
a = num / d
b = y_mean - a * x_mean

在求出a、b之后，可以计算出y的预测值，首先绘制模型直线：

y_hat = a * x + b

plt.scatter(x,y)    # 绘制散点图
plt.plot(x,y_hat,color='r')    # 绘制直线
plt.axis([0,6,0,6])
plt.show()

15724267266219.jpg

然后进行预测：

x_predict = 6
y_predict = a * x_predict + b
print(y_predict)

5.2

2 向量化运算

我们注意到，在计算参数a时：

# a的分子num、分母d
num = 0.0
d = 0.0
for x_i,y_i in zip(x,y):   # zip函数打包成[(x_i,y_i)...]的形式
    num = num + (x_i - x_mean) * (y_i - y_mean)
    d = d + (x_i - x_mean) ** 2
a = num / d

我们发现又这样一个步骤：向量w和向量v，每个向量的对应项，相乘再相加。其实这就是两个向量“点乘”

这样我们就可以使用numpy中的dot运算，非常快速地进行向量化运算。

总的来说：

向量化是非常常用的加速计算的方式，特别适合深度学习等需要训练大数据的领域。

对于 y = wx + b, 若 w, x都是向量，那么，可以用两种方式来计算，第一是for循环：

y = 0
for i in range(n):
    y += w[i]*x[i]
    y += b

另一种方法就是用向量化的方式实现：

y = np.dot(w,x) + b

二者计算速度相差几百倍，测试结果如下：

import numpy as np
import time

a = np.random.rand(1000000)
b = np.random.rand(1000000)

tic = time.time()
c = np.dot(a, b)
toc = time.time()
print("c: %f" % c)
print("vectorized version:" + str(1000*(toc-tic)) + "ms")

c = 0
tic = time.time()
for i in range(1000000):
    c += a[i] * b[i]
toc = time.time()
print("c: %f" % c)
print("for loop:" + str(1000*(toc-tic)) + "ms")

运行结果：

c: 249981.256724
vectorized version:0.998973846436ms
c: 249981.256724
for loop:276.798963547ms

对于独立的样本，用for循环串行计算的效率远远低于向量化后，用矩阵方式并行计算的效率。因此：

只要有其他可能，就不要使用显示for循环。

3 自实现的工程文件

3.1 代码

还记得我们之前的工程文件吗？创建一个SimpleLinearRegression.py，实现自己的工程文件并调用

import numpy as np

class SimpleLinearRegression:
    def __init__(self):
        """模型初始化函数"""
        self.a_ = None
        self.b_ = None

    def fit(self, x_train, y_train):
        """根据训练数据集x_train,y_train训练模型"""
        assert x_train.ndim ==1, \
            "简单线性回归模型仅能够处理一维特征向量"
        assert len(x_train) == len(y_train), \
            "特征向量的长度和标签的长度相同"
        x_mean = np.mean(x_train)
        y_mean = np.mean(y_train)
        num = (x_train - x_mean).dot(y_train - y_mean)  # 分子
        d = (x_train - x_mean).dot(x_train - x_mean)    # 分母
        self.a_ = num / d
        self.b_ = y_mean - self.a_ * x_mean

        return self

    def predict(self, x_predict):
        """给定待预测数据集x_predict，返回表示x_predict的结果向量"""
        assert x_predict.ndim == 1, \
            "简单线性回归模型仅能够处理一维特征向量"
        assert self.a_ is not None and self.b_ is not None, \
            "先训练之后才能预测"
        return np.array([self._predict(x) for x in x_predict])

    def _predict(self, x_single):
        """给定单个待预测数据x_single，返回x_single的预测结果值"""
        return self.a_ * x_single + self.b_

    def __repr__(self):
        """返回一个可以用来表示对象的可打印字符串"""
        return "SimpleLinearRegression()"

3.2 调用

下面我们在jupyter中调用我们自己写的程序：

首先创建一组数据，然后生成SimpleLinearRegression()的对象reg1，然后调用一下

from myAlgorithm.SimpleLinearRegression import SimpleLinearRegression

x = np.array([1.,2.,3.,4.,5.])
y = np.array([1.,3.,2.,3.,5,])
x_predict = np.array([6])
reg = SimpleLinearRegression()
reg.fit(x,y)

输出：SimpleLinearRegression()

reg.predict(x_predict)
reg.a_
reg.a_

输出：
array([5.2])
0.8
0.39999999999999947

y_hat = reg.predict(x)

plt.scatter(x,y)
plt.plot(x,y_hat,color='r')
plt.axis([0,6,0,6])
plt.show()

15724267575477.jpg

4 总结

在本篇文章中，我们实现了简单线性回归算法的代码，并且使用了向量化运算，事实证明，向量化运算能够提高运算效率。

同时我们发现，只要数学公式推导清楚了，实际写代码时没有太多难度的。

那么我们思考一个问题，在之前的kNN算法（分类问题）中，使用分类准确度来评价算法的好坏，那么回归问题中如何评价好坏呢？

公众号二维码.jpeg