A--多元线性回归的python实现

In [1]:

import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt

In [2]:

#创建一个矩阵
a = np.array([[1, 2], [3, 4]])
m = np.mat(a) 
m

[2]:

matrix([[1, 2],
        [3, 4]])

In [4]:

#矩阵运算回顾
# 矩阵转秩
m.T
# 矩阵乘法
m * m 
a * a
# 矩阵⾏列式
np.linalg.det(m)
# 求逆矩阵
m.I
#转换成array
m.A
#降维成一维
m.fattlen

Out[4]:

matrix([[-2. ,  1. ],
        [ 1.5, -0.5]])

 

假设输入数据为DataFrame格式，最后一列为标签值，在此基础上编写线性回归自定义函数(最小二乘）

In [ ]:

#矩阵公式
w=(x.T * x).I * X.T * y 

In [53]:

#根据最小二乘法推导得 w=(x.T * x).I * X.T * y 注：如果(x.T * X)不满足可逆性,那么最小二乘无解，另不满足凸函数，也无解
#又因为特征矩阵在存在多重共线性的情况下，特征矩阵不满足可逆性，所以在做回归之前，需要消除多重共线性
def standRegres(dataSet):
    #把DataFrame转换成array 在转换成matrix，因为DateFrame每一列数据可以不一样，不能直接计算，转换成matirx同时，数据格式也会统一
    xMat = np.mat(dataSet.iloc[:, :-1].values）
    yMat = np.mat(dataSet.iloc[:, -1].values).T
    xTx = xMat.T*xMat
    if np.linalg.det(xTx) == 0:#判断xTx是否是满秩矩阵，若不满秩，则⽆法对其进⾏求逆矩阵的操作
        print("This matrix is singular, cannot do inverse")
        return
    ws = xTx.I * (xMat.T*yMat)
    return ws
#这⾥需要注意的是，当使⽤矩阵分解来求解多元线性回归时，必须添加⼀列全为1的列，⽤于表征线性⽅程截距b。

In [54]:

ex0 = pd.read_table('ex0.txt', header=None)

ex0.head()

Out[54]:

	0	1	2
0	1.0	0.067732	3.176513
1	1.0	0.427810	3.816464
2	1.0	0.995731	4.550095
3	1.0	0.738336	4.256571
4	1.0	0.981083	4.560815

In [55]:

ws = standRegres(ex0)
ws
#返回结果即为各列特征权重，其中数据集第⼀列值均为1，因此返回结果的第⼀个分量表示截距

Out[55]:

matrix([[3.00774324],
        [1.69532264]])

In [56]:

#可视化展示
yhat = ex0.iloc[:, :-1].values * ws
plt.plot(ex0.iloc[:, 1], ex0.iloc[:, 2], 'o')
plt.plot(ex0.iloc[:, 1], yhat)

Out[56]:

[<matplotlib.lines.Line2D at 0x215fd3146a0>]

 

 

模型评价指标残差平⽅和SSE

In [23]:

y = ex0.iloc[:, -1].values
yhat = yhat.flatten()
rss = np.power(yhat - y, 2).sum()
rss

Out[23]:

1.3552490816814904

In [26]:

#将SSE做一个封装
def sseCal(dataSet, regres):#设置参数为 数据集 与 回归方法
    n = dataSet.shape[0] 
    y = dataSet.iloc[:, -1].values
    ws = regres(dataSet)
    yhat = dataSet.iloc[:, :-1].values * ws
    yhat = yhat.reshape([n,])
    rss = np.power(yhat - y, 2).sum()
    return rss
 

In [29]:

sseCal(ex0,standRegres)

Out[29]:

1.3552490816814904

 

模型评价指标 决定系数R_square,决定系数分布在[0, 1]区间内，且越趋近于1，表明拟合程度越好。

In [21]:

sse = sseCal(ex0, standRegres) 
y = ex0.iloc[:, -1].values
sst = np.power(y - y.mean(), 2).sum()
1 - sse / sst

 

Out[21]:

0.9731300889856916

In [31]:

 

#封装R**2
2
def rSquare(dataSet, regres):#设置参数为 数据集 与 回归方法
    sse = sseCal(dataSet, regres) 
    y = dataSet.iloc[:, -1].values
    sst = np.power(y - y.mean(), 2).sum()
    return 1 - sse / sst

 

In [32]:

rSquare(ex0, standRegres)

Out[32]:

0.9731300889856916

 

线性回归的Scikit-Learn实现

In [60]:

from sklearn import linear_model
reg = linear_model.LinearRegression(fit_intercept=True)
reg.fit(ex0.iloc[:, :-1].values, ex0.iloc[:,-1].values)

Out[60]:

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

In [61]:

reg.coef_#返回系数
Out[61]:
array([0.        , 1.69532264])
In [62]:
reg.intercept_#返回截距
Out[62]:
3.0077432426975905

然后计算模型MSE和决定系数

In [63]:

from sklearn.metrics import mean_squared_error, r2_score
yhat = reg.predict(ex0.iloc[:, :-1])
mean_squared_error(y, yhat)

Out[63]:

0.006776245408407454

In [64]:

mean_squared_error(y, yhat)*ex0.shape[0]

Out[64]:

1.3552490816814908

In [65]:

r2_score(y, yhat)

Out[65]:

0.9731300889856916