手写算法-Python代码推广多元线性回归
Python代码实现线性回归一般式的2种方法
1、梯度下降-矩阵形式
上篇文章介绍了一元线性回归,包括Python实现和sklearn实现的实例、对比,以及一些问题点,详情可以看这里:
链接: 手写算法-Python代码实现一元线性回归
里面封装的one_variable_linear()类只适用于一元线性回归,
本篇文章修改代码,推广至多元线性回归,并介绍2种更简洁的方法。
先给大家复习一下矩阵的基本知识:
转置矩阵:
损失函数可表示为:
可以求得:矩阵形式下,偏导的表达式是:
下面附上我的推导证明过程(刚写的):
有了上述表达式,我们修改上次的代码如下:
class linear():
def __init__(self):
pass
#梯度下降法迭代训练模型参数,x为特征数据,y为标签数据,a为学习率,epochs为迭代次数
def fit(self,x,y,a,epochs):
#计算总数据量
m=x.shape[0]
#给x添加偏置项
X = np.concatenate((np.ones((m,1)),x),axis=1)
#计算总特征数
n = X.shape[1]
#初始化W的值,要变成矩阵形式
W=np.mat(np.ones((n,1)))
#X转为矩阵形式
xMat = np.mat(X)
#y转为矩阵形式,这步非常重要,且要是m x 1的维度格式
yMat =np.mat(y.reshape(-1,1))
#循环epochs次
for i in range(epochs):
W=W-a*xMat.T*(xMat*W-yMat)
return W
def predict(self,x,w): #这里的x也要加偏置,训练时x是什么维度的数据,预测也应该保持一样
return np.dot(x,w)
依然用上次的测试数据集,2个代码比较如下:
import numpy as np
import pandas as pd
from sklearn import datasets #sklearn生成数据集都在这里
from matplotlib import pyplot as plt
#生成一个特征的回归数据集
x,y=datasets.make_regression(n_features=1,noise=15,random_state=2020)
plt.scatter(x,y)
plt.show()
class one_variable_linear():
#初始化参数,k为斜率,b为截距,a为学习率,n为迭代次数
def __init__(self,k,b,a,n):
self.k =k
self.b=b
self.a=a
self.n = n
#梯度下降法迭代训练模型参数
def fit(self,x,y):
#计算总数据量
m=len(x)
#循环n次
for i in range(self.n):
b_grad=0
k_grad=0
#计算梯度的总和再求平均
for j in range(m):
b_grad += (1/m)*((self.k*x[j]+self.b)-y[j])
k_grad += (1/m)*((self.k*x[j]+self.b)-y[j])*x[j]
#更新k,b
self.b=self.b-(self.a*b_grad)
self.k=self.k-(self.a*k_grad)
#每迭代10次,就输出一次图像
if i%10==0:
print('迭代{0}'.format(i)+'次')
plt.plot(x,y,'b.')
plt.plot(x,self.k*x+self.b,'r')
plt.show()
self.params= {'k':self.k,'b':self.b}
#输出系数
return self.params
#预测函数
def predict(self,x):
y_pred =self.k * x + self.b
return y_pred
lr=one_variable_linear(k=1,b=1,a=0.1,n=60)
lr.fit(x,y)
旧代码得到的参数如上,
新代码得到的参数如下:
model = linear()
w = model.fit(x,y,a=0.1,epochs=50)
print(w)
我的天,这是什么鬼,怎么和上面的差的这么多(其实是我故意的),明显这个是不正确的,模型完全没有收敛,问题在哪里?
我们细想一下,正常的梯度下降法,前面是带m的,而矩阵形式,我们直接约掉了m,相当于学习率就被放大了m倍,所以这里学习率a应该设置为a/m=0.001,这样a就相等了,迭代次数也相等,冥冥中感觉这次的系数应该是一样的才对。再跑一下代码:
w = model.fit(x,y,a=0.001,epochs=50)
print(w)
哈哈,完全一样,破案了,这里也解释了,之前说的,为什么损失函数前面1/2m这个值,其实对模型的参数没有影响,
但是你的学习率要选择得对,不能可能无法收敛。
现在这个类就是Python线性回归代码的一般式了,设置合理的学习率和迭代次数,就会得到不错的结果。
2、标准方程法
下面来介绍第二种方法,标准方程法。
有了前面的铺垫,这里就很容易理解了,损失函数:
因为这是一个凸函数,因此一定有极小值。根据最小二乘法的原理,我们要对这个损失函数对θ向量求导取0。结果如下式:
这个推导过程中也可以看到,1/2m对最终的系数没有影响,可以直接被约掉。
编写标准方程法代码如下:
class normal():
def __init__(self):
pass
def fit(self,x,y):
m=x.shape[0]
X = np.concatenate((np.ones((m,1)),x),axis=1)
xMat=np.mat(X)
yMat =np.mat(y.reshape(-1,1))
xTx=xMat.T*xMat
#xTx.I为xTx的逆矩阵
ws=xTx.I*xMat.T*yMat
return ws
model =normal()
model.fit(x,y)
求出来的参数为:
这里要注意:XTX的逆矩阵不是什么时候都可以求得出来的,以下情况求不到XTX的逆矩阵:
1、特征数据高度线性相关;
2、n >>m,即特征数量大于样本数量,此时为非满秩矩阵;
sklearn实现对比标准方程法
from sklearn.linear_model import LinearRegression
LR=LinearRegression()
LR.fit(x,y)
LR.intercept_,LR.coef_
和编写的标准方程法得到的参数一模一样,这里回答了之前说过为什么梯度下降法得到的参数和sklearn里面得到的参数不一样的问题,也说明了sklearn中封装的是标准方程法,毕竟真的简单!
下篇介绍非线性回归,当数据表现为非线性时,该怎么处理。
本文地址:https://blog.csdn.net/weixin_44700798/article/details/110452229
下一篇: 这些做妻子的都很厉害,伤不起!