|NO.Z.00010|——————————|BigDataEnd|——|Arithmetic&Machine.v10|——|Machine:监督学习算法.v09|
一、多元线性回归Python实现
Walter Savage Landor:strove with none,for none was worth my strife.Nature I loved and, next to Nature, Art:I warm'd both hands before the fire of life.It sinks, and I am ready to depart ——W.S.Landor
### --- 多元线性回归Python实现
~~~ # 利用矩阵乘法编写回归算法
~~~ 多元线性回归的执行函数编写并不复杂,主要涉及大量的矩阵运算,
~~~ 需要借助Numpy中的矩阵数据格式来完成。首先执行标准化导入:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
~~~ 回顾可能用到的矩阵运算相关知识:
~~~ # 矩阵创建
~~~ numpy中创建矩阵需要使用mat函数,该函数需要输入二维数组。
a = np.array([[1,2],[3,4]])
m = np.mat(a) # NumPy中创建矩阵需要使?mat函数,该函数需要输??维数组
m
~~~ 在Numpy中,矩阵和数组有非常多的相似之处,例如矩阵的索引、切片等常用方法和数组均保持一致。
~~~ # 矩阵运算
~~~ 当然,定义矩阵对象主要是为了能够执行矩阵运算,
~~~ Numpy中有个linalg库专门给矩阵对象提供计算方法,下面简答单回归常用矩阵运算
# 矩阵转置
m.T
# 矩阵乘法
m * m a * a
# 矩阵行列式
np.linalg.det(m)
# 求逆矩阵
m.I
~~~ # 然后编写线性回归函数,同理,我们假设输入数据集为DataFrame,且最后一列为标签值。
def standRegres(dataSet):
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:
print('This matrix is singular,cannot do inverse') #行列式为0,则该矩阵为奇异矩阵,无法求解逆矩阵
return
ws = xTx.I * (xMat.T*yMat)
return ws
~~~ 这里需要提前判断xTx是否是满秩矩阵。若不满秩,则无法对其进行求逆矩阵的操作。
~~~ 定义完函数后,即可测试运行效果,此处我们建立线性随机数进行多元线性回归方程拟合。
~~~ 这里需要注意的是,当使用矩阵分解来求解多元线性回归时,
~~~ 必须添加一列全为1的列,用于表征线性方程截距W0。
rng = np.random.RandomState(1) # 设置随机种子
x = 5*rng.rand(100) # 100个[0,5)的随机数
y = 2*x-5+rng.randn(100) # 真实规律的标签取值
X = pd.DataFrame(x)
Y = pd.DataFrame(y)
ex = pd.DataFrame(np.ones([100,1])) #添加一列权威1的列,表示截距
data = pd.concat([ex,X,Y],axis=1)
~~~ # 数据满足基本的建模要求,然后执行函数运算:
ws = standRegres(data)
ws
~~~ 返回结果即为各列特征权重,其中数据集第一列值均为1,因此返回结果的第一个分量表示截距。
~~~ 然后可用可视化图形展示模型拟合效果∶
yhat = data.iloc[:,:-1].values*ws # 预测标签值
plt.plot(data.iloc[:,1],data.iloc[:,2],'o') #原始数据点
plt.plot(data.iloc[:,1],yhat) # 拟合直线
~~~ # 回顾算法评估指标
~~~ 残差平方和SSE,其计算公式如下所示:
~~~ 其中的m为数据集记录数,yi为实际标签值,y为预测值,则可利用下属表达式进行计算。
y = data.iloc[:,-1].values
yhat = yhat.flatten()
SSE = np.power(yhat-y,2).sum()
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.flatten()
SSE = np.power(yhat-y,2).sum()
return SSE
~~~ # 测试运行:
~~~ 同时,在回归算法中,为了消除数据集规模对残差平方和的影响,往往我们还会计算平均残差MSE:
sseCal(data,standRegres)
~~~ 其中m为数据集样例个数,以及RMSE误差的均方根,为 开平方后所得结果。
~~~ 当然除此之外最常用的还是R=square判定系数,
~~~ 判定系数的计算需要使用之前介绍的组间误差平方和与离差平方和的概念。
~~~ 在回归分析中,SSR表示聚类中类似的组间平方和概念,
~~~ 译为 Sum of squares of the regression,由预测数据与标签均值之间差值的平方和构成。
~~~ 而SST(Total sum of squares)则是实际值和均值之间的差值的平方和:
~~~ 对比之前介绍的聚类分析,此处可对比理解为以点聚点。
~~~ 同样,和轮廓系数类似,最终的判定系数指标也同时结合了 "组内误差" 和 "组间误差" 两个指标。
~~~ # 判定系数R2测度了回归直线对观测数据的拟合程度。
~~~ 若所有观测点都落在直线上,残差平方和为SSE=0,则R2=1,拟合是完全的;
~~~ 如果y的变化与x无关, 完全无助于解释 的变差Y^=Y-, ,则R2=0。
~~~ # 可R2 的取值范围是[0,1]。
~~~ R2越接近1,表明回归平方和占总平方和的比例越大,回归直线与各观测点越接近,
~~~ 用 的变化来解释y值变差( 取值的波动称为变差)的部分就越多,回归直线的拟合程度就越好;
~~~ 反之,R2越接近0,回归直线的拟合程度就越差。
~~~ 接下来,尝试计算上述拟合结果的判定系数。
sse = sseCal(data,standRegres)
y = data.iloc[:,-1].values
sst = np.power(y-y.mean(),2).sum()
1-sse/sst
~~~ 结果为0.91。能够看出最终拟合效果非常好。
~~~ 当然,我们也可编写函数封装计算判断系数的相关操作,同样留?个调整回归函数的接口。
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(data,standRegres)
Walter Savage Landor:strove with none,for none was worth my strife.Nature I loved and, next to Nature, Art:I warm'd both hands before the fire of life.It sinks, and I am ready to depart ——W.S.Landor