线性回归原理以及Python实现

2018.12.9 星期天 阴 biolearn
在统计学中，线性回归是利用线性回归方程对一个或多个自变量和因变量之间的关系进行建模的一种回归分析方法，只有一个自变量的情况称为一元线性回归，大于一个自变量情况的叫做多元线性回归。

线性回归的模型

• 一元线性回归
  

• 考虑两个自变量

• 最终公式可以表达为

• θ 的求解过程

那么如何计算得到参数 θ ？通过线性回归预测得到的预测值和真实值之间都存在一个误差 ε，所以模型可以表示为

由于误差 ε 是独立同分布的，由中心极限定理可得，误差 ε 服从均值为0，方差为某定值 σ² 的高斯分布，可以得到误差 ε 的概率密度函数为

将误差 ε 用真实值减预测值替换可以得到，在给定 x 和 θ 时的关于 y 的概率密度函数

y^（i） 的值彼此之间是独立的，所以将所有的 y 相乘可以得到它的似然函数，因为 x 和 y 是已知的，所以是关于参数 θ 的似然函数为

两边取对数得到

在给定 x 和 θ时，希望预测的 y 值和真实值越接近越好，所以需要求似然函数的最大值，将上述式子的常数项消除，所以变为对下面的式子求取最小值，也就是线性回归中的损失函数

对于由 M 个样本 N 维特征组成的矩阵 X，求解其进行线性回归的参数 θ，上述损失函数可以变成

其中 Xθ 理解为预测值，y 表示真实值，为了求解上述式子的最小值，根据最小值必定为函数的驻点，可以对上述式子进行求导

驻点的导数为0，于是得到了参数 θ 的解析式为

以上是参数 θ 的解析式的推导和求解过程，是利用最小二乘法的原理完成的，但是由于求解矩阵的逆很耗费时间，一般是通过梯度下降算法进行计算的。

Python实现线性回归

根据上述推到的求解参数的公式计算参数，构建模型并进行预测，将结果和利用Scikit-Learn中线性回归得到的结果进行比较，查看参数和预测结果是否一致

• 数据描述

数据来自出版书籍《An Introduction to Statistical Learning with Applications in R》 (Springer, 2013)，作者Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani。共200 条数据， 每条数据4 个属性。该数据可以从这个链接直接下载得到 http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv

• 数据集信息，数据共4列200行，每一行为一个特定的商品，前三列为输入特征也就是x，最后一列为输出特征也就是y

• python

#!/usr/bin/python
# -*- coding:utf-8 -*-
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats
from sklearn import linear_model

if __name__ == "__main__":
    data = pd.read_csv('./Advertising.csv')    # TV,Radio,Newspaper,Sales
    data['intercept'] = 1
    x = data[['intercept','TV', 'Radio', 'Newspaper']]
    y = data['Sales']
    # 前150行的数据作为train，后50行作为test
    train_x = np.array(x.loc[1:150,])
    test_x = np.array(x.loc[151:,])
    train_y = np.array(y.loc[1:150,])
    test_y = np.array(y.loc[151:,])
    # beta = (X^TX)^(-1)X^Ty，计算参数
    Xt = np.transpose(train_x)
    XtX = np.dot(Xt,train_x)
    Xty = np.dot(Xt,train_y)
    beta = np.linalg.solve(XtX,Xty)
    print(beta)
    # 对后50行的数据进行预测
    pred=[]
    for data, actual in zip(test_x, test_y):
        test = np.transpose(data)
        prediction = np.dot(test, beta)
        pred.append(prediction)
        print('prediction = ' + str(prediction) + ' actual = ' + str(actual))
    #plot
    pcc = stats.pearsonr(test_y, pred)[0]
    print('Pearson Correlation Coefficient = ' + str(pcc))

    mpl.rcParams['font.sans-serif'] = ['simHei']
    mpl.rcParams['axes.unicode_minus'] = False
    plt.figure(figsize=(14, 6),facecolor='w')
    plt.subplot(121)
    t = np.arange(len(test_x))
    plt.plot(t, test_y, 'r-', linewidth=2, label='真实数据')
    plt.plot(t, pred, 'g-', linewidth=2, label='预测数据')
    plt.title('线性回归预测销量', fontsize=18)
    plt.text(10, 20, 'PCC=%0.2f' % pcc, fontsize=15)
    plt.legend(loc='upper left')
    plt.grid(b=True, ls=':')

    ## 利用sklearn中的linear_model
    model = linear_model.LinearRegression()
    model.fit(train_x[:,1:],train_y)
    print('利用sklearn中的linear_model计算的参数：%s\t%s' % (model.coef_, model.intercept_))
    test_y_pred = model.predict(test_x[:,1:])
    pcc_sk = stats.pearsonr(test_y,test_y_pred)[0]
    print('Pearson Correlation Coefficient(sklearn) = ' + str(pcc_sk))
    plt.subplot(122)
    plt.plot(t, test_y, 'r-', linewidth=2, label='真实数据')
    plt.plot(t, pred, 'go-', linewidth=2, label='预测数据')
    plt.title('线性回归预测销量(sklearn)', fontsize=18)
    plt.text(10, 20, 'PCC=%0.2f' % pcc_sk, fontsize=15)
    plt.legend(loc='upper left')
    plt.grid(b=True, ls=':')
    plt.show()

• 结果输出，可以看到通过自行编写得到的参数以及预测结果和通过使用sklearn得到的结果一模一样

[ 3.07875053e+00  4.65616125e-02  1.80900459e-01 -2.55988893e-03]   #截距以及三个特征的权重
Pearson Correlation Coefficient = 0.9505078535320115
利用sklearn中的linear_model计算的参数：[ 0.04656161  0.18090046 -0.00255989]  3.0787505315686783  # 三个特征的权重以及截距
Pearson Correlation Coefficient(sklearn) = 0.9505078535320114