线性回归之最小二乘法
线性回归
Linear Regression——线性回归
是机器学习中有监督机器学习下的一种简单的回归算法。
分为一元线性回归(简单线性回归)和多元线性回归,其中一元线性回归是多元线性回归的一种特殊情况,我们主要讨论多元线性回归如果因变量和自变量之间的关系满足线性关系(自变量的最高幂为一次),那么我们可以用线性回归模型来拟合因变量与自变量之间的关系.
简单线性回归的公式如下:
y^=ax+b hat y=ax+by^?=ax+b
多元线性回归的公式如下:
y^=θTx hat y= theta^T x y^?=θTx
上式中的θthetaθ为系数矩阵,x为单个多元样本.
由训练集中的样本数据来求得系数矩阵,求解的结果就是线性回归模型,预测样本带入x就能获得预测值y^hat yy^?,求解系数矩阵的具体公式接下来会推导.
推导过程
推导总似然函数
假设线性回归公式为y^=θxhat y= theta xy^?=θx.
真实值y与预测值y^hat yy^?之间必然有误差?=y^?yepsilon=hat
y-y?=y^?y,按照中心极限定理(见知识储备),我们可以假定?epsilon?服从正态分布,正态分布的概率密度公式为:
ρ(x)=1σ2πe?(x?μ)22σ2rho (x)=frac {1}{sigmasqrt{2pi}}e^{-frac{(x-mu)^2}{2sigma^2}}ρ(x)=σ2π?1?e?2σ2(x?μ)2?
为了模型的准确性,我们希望?epsilon?的值越小越好,所以正态分布的期望μmuμ为0.概率函数需要由概率密度函数求积分,计算太复杂,但是概率函数和概率密度函数呈正相关,当概率密度函数求得最大值时概率函数也在此时能得到最大值,因此之后会用概率密度函数代替概率函数做计算.
我们就得到了单个样本的误差似然函数(μ=0,σmu=0,sigmaμ=0,σ为某个定值):
ρ(?)=1σ2πe?(?0)22σ2rho (epsilon)=frac {1}{sigmasqrt{2pi}}e^{-frac{(epsilon-0)^2}{2sigma^2}}ρ(?)=σ2π?1?e?2σ2(?0)2?
而一组样本的误差总似然函数即为:
Lθ(?1,?,?m)=f(?1,?,?m∣μ,σ2)L_theta(epsilon_1,cdots,e psilon_m)=f(epsilon_1,cdots,epsilon_m|mu,sigma^2)Lθ?(?1?,? ,?m?)=f(?1?,?,?m?∣μ,σ2)
因为我们假定了?epsilon?服从正态分布,也就是说样本之间互相独立,所以我们可以把上式写成连乘的形式:
f(?1,?,?m∣μ,σ2)=f(?1∣μ,σ2)?f(?m∣μ,σ2)f(epsilon_
1,cdots,epsilon_m|mu,sigma^2)=f(epsilon_1|mu,sigma^2)*cdots *f(epsilon_m|mu,sigma^2)
f(?1?,?,?m?∣μ,σ2)=f(?1?∣μ,σ2)?f(?m?∣μ,σ2) Lθ(?1,?,?m)=∏i=1mf(?i∣μ,σ2)=∏i=1m1σ2πe?(?i?0)22σ2L_theta(epsilon_1,cdots,epsilon_m)=prod^m_{i=1}f(epsilon _i|mu,sigma^2)=prod^m_{i=1}frac
{1}{sigmasqrt{2pi}}e^{-frac{(epsilon_i-0)^2}{2sigma^2}}Lθ? (?1?,?,?m?)=i=1∏m?f(?i?∣μ,σ2)=i=1∏m?σ2π?1?e?2σ2(?i?
0)2?
在线性回归中,误差函数可以写为如下形式:
?i=∣yi?y^i∣=∣yi?θTxi∣epsilon_i=|y_i-hat
y_i|=|y_i-theta^Tx_i|?i?=∣yi?y^?i?∣=∣yi?θTxi?∣最后可以得到在正态分布假设下的总似然估计函数如下:
Lθ(?1,?,?m)=∏i=1m1σ2πe?(?i?0)22σ2=∏i=1m1σ2πe?(yi ?θTxi)22σ2L_theta(epsilon_1,cdots,epsilon_m)=prod^m_{i=1} frac
{1}{sigmasqrt{2pi}}e^{-frac{(epsilon_i-0)^2}{2sigma^2}}=pro d^m_{i=1}frac
{1}{sigmasqrt{2pi}}e^{-frac{(y_i-theta^Tx_i)^2}{2sigma^2}}L θ?(?1?,?,?m?)=i=1∏m?σ2π?1?e?2σ2(?i?0)2?=i=1∏m?σ2π?1 ?e?2σ2(yi?θTxi?)2?
推导损失函数
按照最大总似然的数学思想(见知识储备),我们可以试着去求总似然的最大值.
遇到连乘符号的时候,一般思路是对两边做对数运算(见知识储备),获得对数总似然函数:
l(θ)=loge(Lθ(?1,?,?m))=loge(∏i=1m1σ2πe?(yi?θTxi)22σ2)l(theta)=log_e(L_theta(epsilon_1,cdots,epsilon_m))=log_ e(prod^m_{i=1}frac
{1}{sigmasqrt{2pi}}e^{-frac{(y_i-theta^Tx_i)^2}{2sigma^2}}) l(θ)=loge?(Lθ?(?1?,?,?m?))=loge?(i=1∏m?σ2π?1?e?2σ2(yi ?θTxi?)2?)
l(θ)=loge(∏i=1m1σ2πe?(yi?θTxi)22σ2)=∑i=1mloge1σ2πexp(?(yi?θTxi)22σ2)=mloge1σ2π?12σ2∑i=1m(yi?θTxi)2l (theta) = log_e(prod^m_{i=1}frac {1}{sigmasqrt{2pi}}e^{-frac{(y_i-theta^Tx_i)^2}{2sigma^2}}) = sum_{i=1}^mlog_efrac {1}{sigmasqrt{2pi}}exp({-frac{(y_i-theta^Tx_i)^2}{2sigma^2} })=mlog_efrac{1}{sigmasqrt{2pi}}-frac{1}{2sigma^2}sum^m_{i= 1}(y^i-theta^Tx^i)^2l(θ)=loge?(i=1∏m?σ2π?1?e?2σ2(yi?θTxi?)2?)=i=1∑m?loge?σ2π?1?exp(?2σ2(yi?θTxi?)2?)=mloge?σ2π?1?2σ21?i=1∑m?(yi?θTxi)2
前部分是一个常数,后部分越小那么总似然值越大,后部分则称之为损失函数,则有损失函数的公式J(θ)J(theta)J(θ):
J(θ)=12∑i=1m(yi?θTxi)2=12∑i=1m(yi?hθ(xi))2=12∑i=1m (hθ(xi)?yi)2
J(theta)=frac{1}{2}sum^m_{i=1}(y^i-theta^Tx^i)^2=frac{1}{2} sum^m_{i=1}(y^i-h_theta(x^i))^2=frac{1}{2}sum^m_{i=1}(h_the ta(x^i)-y^i)^2
J(θ)=21?i=1∑m?(yi?θTxi)2=21?i=1∑m?(yi?hθ?(xi))2=21?i=1∑m?(hθ?(xi)?yi)2
解析方法求解线性回归
要求的总似然最大,需要使得损失函数最小,我们可以对损失函数求导.
首先对损失函数做进一步推导:
J(θ)=12∑i=1m(hθ(xi)?yi)2=12(Xθ?y)T(Xθ?y)J(theta)=fr ac{1}{2}sum^m_{i=1}(h_theta(x^i)-y^i)^2=frac{1}{2}(Xtheta-y )^T(Xtheta-y)J(θ)=21?i=1∑m?(hθ?(xi)?yi)2=21?(Xθ?y)T(Xθ?y)
注意上式中的X是一组样本形成的样本矩阵,θthetaθ是系数向量,y也是样本真实值形成的矩阵,这一步转换不能理解的话可以试着把
12(Xθ?y)T(Xθ?y)frac{1}{2}(Xtheta-y)^T(Xtheta-y)21?(Xθ?y) T(Xθ?y)带入值展开试试.
J(θ)=12∑i=1m(hθ(xi)?yi)2=12(Xθ?y)T(Xθ?y)=12((Xθ)T? yT)(Xθ?y)=12(θTXT?yT)(Xθ?y)=12(θTXTXθ?yTXθ?θTXTy+yTy
)J(theta)=frac{1}{2}sum^m_{i=1}(h_theta(x^i)-y^i)^2=frac{1} {2}(Xtheta-y)^T(Xtheta-y)=frac{1}{2}((Xtheta)^T-y^T)(Xtheta -y)=frac{1}{2}(theta^TX^T-y^T)(Xtheta-y)=frac{1}{2}(theta^T X^TXtheta-y^TXtheta-theta^TX^Ty+y^Ty)J(θ)=21?i=1∑m?(hθ?( xi)?yi)2=21?(Xθ?y)T(Xθ?y)=21?((Xθ)T?yT)(Xθ?y)=21?(θTXT ?yT)(Xθ?y)=21?(θTXTXθ?yTXθ?θTXTy+yTy)
根据黑塞矩阵可以判断出J(θ)J(theta)J(θ)是凸函数,即J(θ)J(theta)J(θ)的对θthetaθ的导数为零时可以求得J(θ)J(theta)J(θ)的最小值.
?J(θ)?θ=12(2XTXθ?(yTX)T?XTy)=12(2XTXθ?XTy?XTy)=XTXθ?XTyfrac{partial
J(theta)}{partialtheta}=frac{1}{2}(2X^TXtheta-(y^TX)^T-X^Ty )=frac{1}{2}(2X^TXtheta-X^Ty-X^Ty)=X^TXtheta-X^Ty?θ?J(θ)? =21?(2XTXθ?(yTX)T?XTy)=21?(2XTXθ?XTy?XTy)=XTXθ?XTy 当上式等于零时可以求得损失函数最小时对应的θthetaθ,即我们最终想要获得的系数矩阵:
XTXθ?XTy=0XTXθ=XTy((XTX)?1XTX)θ=(XTX)?1XTyEθ=(XTX)?1 XTyθ=(XTX)?1XTyX^TXtheta-X^Ty=0X^TXtheta=X^Ty
((X^TX)^{-1}X^TX)theta=(X^TX)^{-1}X^Ty
Etheta=(X^TX)^{-1}X^Tytheta=(X^TX)^{-1}X^TyXTXθ?XTy=0XT Xθ=XTy((XTX)?1XTX)θ=(XTX)?1XTyEθ=(XTX)?1XTyθ=(XTX)?1XTy (顺便附上一元线性回归的系数解析解公式:θ=∑i=1m(xi?x
￣)(yi?y￣)∑i=1m(xi?x ￣)2theta=frac{sum^m_{i=1}(x_i-overline{x})(y_i-overline{y} )}{sum^m_{i=1}(x_i-overline{x})^2}θ=∑i=1m?(xi?x)2∑i=1m?( xi?x)(yi?y?)?)
简单实现
import numpy as np
import matplotlib.pyplot as plt
# 随机创建训练集,X中有一列全为'1'作为截距项
X = 2 * np.random.rand(100, 1)
y = 5 + 4 * X + np.random.randn(100, 1)
X = np.c_[np.ones((100,1)),X]
# 按上面获得的解析解来求得系数矩阵theta
theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
# 打印结果
print(theta)
# 测试部分
X_test = np.array([[0],
X_test = np.c_[(np.ones((2, 1))), X_test]
print(X_test)
y_predict = X_test.dot(theta)
print(y_predict)
plt.plot(X_test[:,-1], y_predict, 'r-')
plt.axis([0, 2, 0, 15])
plt.show()
sklearn实现
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression X = 2 * np.random.rand(100, 1)
y = 5 + 4 * X + np.random.randn(100, 1)
X = np.c_[np.ones((100,1)),X]
# 新建线性回归模型
model = LinearRegression(fit_intercept=False)
# 代入训练集数据做训练
model.fit(X,y)
# 打印训练结果
print(model.intercept_,model.coef_)
X_test = np.array([[0],
X_test = np.c_[(np.ones((2, 1))), X_test]
print(X_test)
y_predict =model.predict(X_test)
print(y_predict)
plt.plot(X_test[:,-1], y_predict, 'r-')
plt.axis([0, 2, 0, 15])
plt.show()
使用解析解的公式来求得地模型是最准确的.
计算量非常大,这会使得求解耗时极多,因此我们一般用的都是梯度下降法求解.
知识储备
距离公式
机器学习中常见的距离公式 - WingPig - 博客园
中心极限定理
是讨论随机变量序列部分和分布渐近于正态分布的一类定理。这组定理是数理统计学和误差分析的理论基础，指出了大量随机变量累积分布函数逐点收敛到正态分布的积累分布函数的条件。线性回归中的误差?epsilon?就可以看作是一个随机变量.
在自然界与生产中，一些现象受到许多相互独立的随机因素的影响，如果每个因素所产生的影响都很微小时，总的影响可以看作是服从正态分布的,中心极限定理从数学上证明了这一现象.
按照中心极限定理,我们可以认为当随机变量?epsilon?足够多的时候,是服从正态分布(即高斯分布).正态分布的概率密度函数公式为:
f(x)=1σ2πe?(x?μ)22σ2f(x)=frac
{1}{sigmasqrt{2pi}}e^{-frac{(x-mu)^2}{2sigma^2}}f(x)=σ2π?
1?e?2σ2(x?μ)2?
机器学习中我们常常假设误差符合均值为0，方差为定值的正态分布！！
最大似然值估计
是一种统计方法
求一个样本集的相关概率密度函数的参数,即估计一个数据集的概率密度函数.在线性回归算法中,要做的就是求正态分布函数中的μmuμ,这里用来估计正态分布中的方差参数值.
我们要求得的是似然值最大时的参数.
假设每一个误差?epsilon?都服从正态分布,将所有的?epsilon?叠加起来可以获得总似然,即将每一个?epsilon?的概率函数相乘.
但是要求正态分布的概率函数需要对概率密度函数做积分,为了便于运算这里可以直接把正态分布的概率密度函数相乘,找到概率密度相乘最大时相当于找到了概率相乘最大的时候.
矩阵求导
常用的矩阵求导公式有:
因为我们要根据样例数据去预测未知结果，所以我们的样例数据格式为：
#xA0;#xA0;#xA0;#xA0;#xA0;#xA0;" role="presentation" style="position: relative;"> ~~~~~~for j = 1 to n{ end{aligned}?w?(yTXw)?=(yTX)T=XTy，分母布局?
我们的目标便是选择出可以使得建模误差的平方和能够最小的模
型参数。即使得损失函数J(θ)J(theta)J(θ)最小，也就是说，我们通过这个方式，能够得到θ0,θ1theta_0, theta_1θ0?,θ1?的合适大小
from sklearn.linear_model import LinearRegression
#可使用loadtxt()函数读取txt文档，delimiter=','代表以逗号分隔开
注: 此时的整合中的x = (x0 x1 x2 . xn), x0是新加入的值，为了保证θ0起作用，则x0 = 1。为什么要进行加入x0这步操作呢？其实就是为了化简表达式, 使的可以使用矩阵来进行表示。
令 v?=Xu?vec{v} = Xvec{u}v=Xu :
(XTX)?1(XTX)w=(XTX)?1XTy(X^TX)^{-1}(X^TX)w=(X^TX)^{-1}X^ Ty(XTX)?1(XTX)w=(XTX)?1XTy
我们希望求出的W是最接近线性方程的解的，最接近我们定义为残差平方和最小，残差的公式和残差平方和的公式如下：