什么是多元线性回归?
在线性回归分析中,如果有两个或两个以上的自变量,就称为多元线性回归(multivariable linear regression)。如果我们要预测一套房子的价格,影响房子价格的因素可能包括:面积、卧室数量、层数以及房龄,我们用x1、x2、x3、x4来代表这4个特征(features)。在这里,我们假设房子的价格和这4个特征值都是线性相关的,并且用hθ来表示预测房子的价格,那么有:
hθ(x)=θ0+θ1x1+θ2x2+θ3x3+θ4x4
其中θ为线性回归中的参数。我们把x、θ分别表示为列向量的形式:
θ=(θ0,θ1,θ2,θ3,θ4)T
x=(1,x1,x2,x3,x4)T
那么我们可以发现
hθ(x)=θT∙x
如何实现多元线性回归?
对于拥有n个特征值、m个数据样本的数据,我们可以用一个m*(n+1)矩阵的形式表示出来,其中矩阵的每一行(即x(i))为一个数据样本,每一列代表一个特征值,xj(i)代表第i个数据样本中第j个特征值(x0(i)=1)
X=⎛⎝⎜⎜⎜⎜⎜⎜⎜111⋮1x(1)1x(2)1x(3)1⋮x(m)1x(1)2x(2)2x(3)2⋮x(m)2⋯⋯⋯⋯⋯x(1)nx(2)nx(3)n⋮x(m)n⎞⎠⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜x(1)Tx(2)Tx(3)T⋮x(m)T⎞⎠⎟⎟⎟⎟⎟⎟⎟
样本的实际值(对于房价预测来说就是房价)用y来表示:
y=(y1,y2,…,ym)T
线性回归参数θ为:
θ=(θ0,θ1,…,θn)T
接下来我们定义代价函数Jθ(cost function)为:
Jθ(X)=12m∑i=1m(hθ(x(i))−y(i))2=12m∑i=1m(θT∙x(i)−y(i))2=12m∑i=1m(x(i)T∙θ−y(i))2
我们还可以把Jθ利用矩阵表达,获得更为简洁的形式
Jθ(X)=12m(X∙θ−y)T∙(X∙θ−y)
现在,我们只需要让代价函数Jθ最小,就能得到最优的θ参数。那么,要怎样才能使Jθ最小呢?有两个办法,一个是梯度下降法(gradient descent),一个是标准方程法(norm equation)。
Jθ在样本数据X确定时,是关于θ的一个函数,而θ是一个(n+1)维的列向量,也就是说Jθ其实是一个(n+1)元函数。学过微积分的人应该知道,要求一个多元函数的极值,可以令多元函数对每个未知元的偏导数为0,从而求出未知元,再代入函数中,便可求出函数的极值。同样地,我们现在要求Jθ的极值,那么就要先使Jθ关于θ的偏导数,或者说关于θ0、θ1、θ2、…、θn的偏导数都为0。
Jθ对θj的偏导数为:
δJδθj=1m∑i=1m(x(i)T∙θ−y(i))x(i)j=1m(x(1)jx(2)j⋯x(m)j)∙(Xθ−y)
则Jθ对θ的偏导数为:
δJδθ=⎛⎝⎜⎜⎜⎜⎜⎜δJδθ0δJδθ1⋮δJδθn⎞⎠⎟⎟⎟⎟⎟⎟=1m⎛⎝⎜⎜⎜⎜⎜1x(1)1⋮x(1)n1x(2)1⋮x(2)n⋯⋯1x(m)1⋮x(m)n⎞⎠⎟⎟⎟⎟⎟(Xθ−y)=1mXT(Xθ−y)
梯度下降法
梯度下降法的思想是一点点、逐步靠近极小值,通过多次迭代,不断更新θ的值,让代价函数Jθ收敛于极小值。迭代的公式为:
上述理论参考:https://www.cnblogs.com/magic-girl/p/mutivariable-linear-regression.html
代码实现:
准备数据:
可以将以上数据保存在Delivery.csv中。
我们想要根据以上数据预测103, 7 所对应的y值。
from numpy import genfromtxt
import numpy as np
from numpy.linalg import pinv
dataPath = r"Delivery.csv"
# 加上r就会将后面当做一个完整的字符串,而不会认为里面有什么转义之类的。
deliveryData = genfromtxt(dataPath, delimiter=',')
print("data")
print(deliveryData)
x = deliveryData[:, :-1]
y = deliveryData[:, -1]
print("x: ")
print(x)
print("y: ")
print(y)
#一维numpy数组转置方法
y = y.reshape(len(y),-1)
print(y)
# 标准方程法(Normal Equation)
# theta = pinv(x.T.dot(x)).dot(x.T).dot(y)
# 梯度下降法(Gradient Descent)
#δJδθ=1/mXT(Xθ−y)
def partial_derivative(X, theta, y):
print(X.shape, theta.shape, y.shape)
print(X.T.dot(X.dot(theta) - y))
derivative = X.T.dot(X.dot(theta) - y) / X.shape[0]
print(derivative)
return derivative
def gradient_descent(X, y, alpha=0.0001):
print(X.shape[1])
theta = np.ones(shape=X.shape[1], dtype=float)
theta = theta.reshape(len(theta), -1)
partial_derivative_of_J = partial_derivative(X, theta, y)
while any(abs(partial_derivative_of_J) > 1e-5):
theta = theta - alpha * partial_derivative_of_J
print("theta.shape", theta.shape)
partial_derivative_of_J = partial_derivative(X, theta, y)
return theta
theta = gradient_descent(x, y)
print(theta)
# print(theta.shape)
# 预测一个例子
xPred = np.array([103, 7])
# print(xPred.shape)
print(xPred)
# hθ(x)=θT∙x
yPred = theta.T.dot(xPred)
print("predicted y:")
print(yPred)
#[ 11.32750806]
运行后的结果为11.32750