前言
一元线性回归只是让我们好理解回归这个概念。在实际中特征很少说只有一个。所有就有多元线性回归。
概念
与一元线性回归一样,多元线性回归自然是一个回归问题。只不过一元线性回归的方程可以写成。
y=ax+b
多元线性回归是
y=a1x1+a2X2+a3X3+...+anXn+b
相当于我们高中学的一元一次方程,变成了n元一次方程。因为y还是那个y。只是特征增多了。
模型的训练
多元线性回归的训练方式与一元线性回归方式基本相同。
1、设置一个损失函数。
2、想办法让润湿函数最小。
其中
每一个i对应一行的数据。
我们只需要找到合适的。
参数即可。
方便计算
为了方便算,我们把b换一下。
那么可以组成一个待解参数矩阵。
然而我们的输入X对应每个特征也可以组成一个矩阵。
于是我们很快发现 X的每一行与O不同维。O的维度是1*n+1 (从0开始数起的)
X的一行的维度是是1Xn
那么我们只有在X矩阵加一列。那么加什么不会改变O1的值呢,那就是1了。
那么我们X矩阵就变成了
大家都知道矩阵的乘法是行*列然后每个元素想加。这里的yi就等于Xi 与 O的点乘,因为O是一个行向量,所以这里要转置成列向量。
(备注:numpy.array 中的dot方法会自动转换行、列向量,但是这里说原理要说清楚)。
那么可以得到
其中Xb就是X加一列的结果。O是一个列向量。所得到的y的预测值也是一个向量。
损失函数
我们再看损失函数,长这样。
其中的y是一个1 X m的向量,其中每一元素代表是一个标签。y预测也是一个1 * m的向量。这个公式计算每一个样本的误差然后求和。那么我们可以直接转换成矩阵的计算。
我们可以改为
我们只需要找到 合适的O就可以了。
找参数
我们只需要找到合适的cita就可以了。
正规方程解
正规方程解其实就是跟简单线性回归一样求偏导数然后求极值点。这里多元的那就是求偏导数然后求极值点,但是多元的好麻烦。这里直接给出公式比较合适。
这个知道仔细推导其实意义并不大,因为我们在实际算法中不会去使用这个。我们都知道矩阵的乘法时间复杂度是O(n3)有很多人做了优化也没有降到O(n2),如果你能降低0.1个点就不得了。所以说这是一个世界性难题。大家都知道我们机器学习其实计算量是很大的,那么怎么能让我们的计算变的快一点只有改变我们的策略。那就是梯度下降。
使用求正规方程解寻找参数代码
def fit_normal(self,X_train,y_train):
assert X_train.shape[0] == y_train.shape[0],"你要一一对应吧"
X_b = np.hstack([np.ones((len(X_train),1)),X_train])
self._theta=np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train) #就是那个公式
self.interception_=self._theta[0] #第0项就是那个截距b
self.coef_=self._theta[1:] #后面都是系数
梯度下降
在简单线性回归的时候说了梯度下降的思路。其实就是一种求极值的一种思路,就是在模型上随机选一个点,然后向着梯度方向以一定步长挪动。最终能到达极值点。
在简单线性回归中就是x沿着x的导数这个梯度进行下降。
可以把loss函数定义一下
那么两个参数的变化是
同样的多元线性回归是一样的。其中参数都是cta
可以定义梯度
那么每次更新就是
那么关键就是求梯度这个矩阵。
可以看一下。
注意这里的X 是在第一列加了一个全1项的。
那么梯度矩阵可以写为
我们会发现一个问题就是m越大,梯度就越大,梯度怎么会跟m(样本的数量有关系呢,这显然不合理)。所以损失函数应该 加一个1/m这个常数来抵消这个影响。
那么梯度就变成了。
看代码吧
def J(theta,X_b,y):
"""
定义损失函数
"""
try:
return np.sum((y- X_b.dot(theta))**2)/len(X_b) #=损失函数
except:
return float('inf')
def dJ(theta,X_b,y):
"""
求偏导,计算梯度
"""
res = np.empty(len(theta))
res[0]=np.sum(X_b.dot(theta)-y)
for i in range(1,len(theta)):
res[i]=(X_b.dot(theta)-y).dot(X_b[:,i])
return res * 2/len(X_b)
def gradient_descent(X_b,y,initial_theta,eta,n_iters=1e4,epsilon=1e-8):
theta=initial_theta
i_iter=0
while i_iter < n_iters:#防止梯度不停的跳动,过了一定轮数就停止。
gradient = dJ(theta,X_b,y)
last_theta= theta
theta=theta - eta * gradient
if(abs(J(theta,X_b,y)-J(last_theta,X_b,y))<epsilon):# 不一定非要是0,接近0就可以了。
break
i_iter +=1
return theta