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]])

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

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

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

#矩阵公式
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

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

Out[55]:

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

#可视化展示
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

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

1.3552490816814904

#将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，表明拟合程度越好。

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

 

0.9731300889856916

 

#封装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

 

rSquare(ex0, standRegres)

Out[32]:

0.9731300889856916

线性回归的Scikit-Learn实现

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)

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

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

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

0.006776245408407454

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

Out[64]:

1.3552490816814908

r2_score(y, yhat)

Out[65]:

0.9731300889856916