线性回归-2-简单线性回归的代码实现

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

上篇文章简单介绍了简单线性回归，并从数学角度加以理解，本篇文章使用python实现一个简单的线性回归小例子。

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

1.1、数据源引入

首先引入一个数据源,该数据源中保存了工作年限与薪水的若干条相关信息

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

income = pd.read_csv(r'./file/Salary_Data.csv')
income.head()

# 绘制散点图
sns.lmplot(x = 'YearsExperience'
           , y = 'Salary'
           , data = income
           , ci = None
           , fit_reg=True)
# 显示图形
plt.show()

output_6_1.png

这里使用seaborn可以直接得到拟合的直线，这条直线就是我们需要预测的。

1.2、思路梳理

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

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

在本篇文章中，我们直接假设模型是一条直线，模型是： $\sum_{i=1}^{m}\left(y^{(i)}-\hat{y}^{(i)}\right)^{2}$

根据最小二乘法推导求出a、b的表达式：
$a=\frac{\sum_{i=1}^{m}\left(x^{i}-\bar{x}\right)\left(y^{i}-\bar{y}\right)}{\sum_{i=1}^{m}\left(x^{(i)}-\bar{x}\right)^{2}} \quad b=\bar{y}-a \bar{x}$

1.3、编码实现

首先我们计算x、y的均值,通过均值计算a、b的值

x_mean = round(income.YearsExperience.mean(),2)
y_mean = income.Salary.mean()
x_mean
y_mean
a_upper = 0.0
a_lower = 0.0

# for i,k in enumerate(income.YearsExperience):
#     a_upper += (k - x_mean)*(income.Salary[i] - y_mean)
#     a_lower += (k - x_mean)**2

# 或者写成一下形式，同时思考一下，在大量数据的情况下有何优化措施？
for x_i,y_i in zip(income.YearsExperience, income.Salary):   # zip函数打包成[(x_i,y_i)...]的形式
    a_upper = a_upper + (x_i - x_mean) * (y_i - y_mean)
    a_lower = a_lower + (x_i - x_mean) ** 2
    
a = a_upper/a_lower
b = y_mean - a * x_mean
a
b

5.31


76003.0


9449.94883432168


25823.771689751884

我们使用拟合好的模型预测样本中的y值并直接可视化，看看与seaborn中的是否相似

import matplotlib.pyplot as plt
%matplotlib inline

y_pre = a*income.YearsExperience + b
plt.figure(figsize=(4^2,4^2))
plt.scatter(income.YearsExperience, income.Salary)    # 绘制散点图
plt.plot(income.YearsExperience, y_pre,color='r')    # 绘制直线

2、numpy向量化运算

2.1、向量化运算简介与使用

我们可以使用numpy中的dot运算，非常快速地进行向量化运算。下面来介绍一下dot的使用：

dot()返回的是两个数组的点积(dot product)

如果处理的是一维数组，则得到的是两数组的內积

import numpy as np
d = np.arange(0,9)
# array([0, 1, 2, 3, 4, 5, 6, 7, 8])
e = d[::-1]
# array([8, 7, 6, 5, 4, 3, 2, 1, 0])
np.dot(d,e)

如果是二维数组（矩阵）之间的运算，则得到的是矩阵积（mastrix product）。

a = np.arange(1,5).reshape(2,2)

# array([[1, 2],
#      [3, 4]])

b = np.arange(5,9).reshape(2,2)
# array([[5, 6],
#       [7, 8]])
np.dot(a,b)

array([[19, 22],
       [43, 50]])

所得到的数组中的每个元素为：第一个矩阵中与该元素行号相同的元素与第二个矩阵与该元素列号相同的元素，两两相乘后再求和。下面使用一张图片理解这句话。

点积.png

此时我们再来看我们之前写的代码：

for x_i,y_i in zip(income.YearsExperience, income.Salary):   # zip函数打包成[(x_i,y_i)...]的形式
    a_upper = a_upper + (x_i - x_mean) * (y_i - y_mean)
    a_lower = a_lower + (x_i - x_mean) ** 2

a_upper = a_upper + (x_i - x_mean) * (y_i - y_mean)

实际上就是两数组的内积，
可以写成

(income.YearsExperience - x_mean).dot(income.Salary - y_mean)

2.2、向量化运算的优点

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

以本文中的数据为例，比较向量化运算与普通计算速度的差异：

import time

a = income.YearsExperience
b = income.Salary
for i in range(500):
    a = a.append(income.YearsExperience)
    b = b.append(income.Salary)
a = a.reset_index()
b = b.reset_index()
tic = time.time()
c = np.dot(a.YearsExperience, b.Salary)
toc = time.time()
# print(len(a))
print("c: %f" % c)
print("vectorized version:" + str(1000*(toc-tic)) + "ms")

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

c: 7175302461.000004
vectorized version:0.0ms
c: 7175302460.999892
for loop:752.0430088043213ms

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

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

3、简单的代码封装与调用

3.1、封装

%%writefile 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()"

Overwriting SimpleLinearRegression.py

3.2、调用

from SimpleLinearRegression import SimpleLinearRegression

x = income.YearsExperience
y = income.Salary

reg = SimpleLinearRegression()
reg.fit(x,y)


reg.a_
reg.b_

9449.962321455077


25792.20019866869

x_predict = np.array([6])

reg.predict(x_predict)

array([82491.9741274])