您现在的位置是:网站首页> 内容页

多元线性回归公式推导及R语言实现

  • www.999dabao.com
  • 2019-07-19
  • 32人已阅读
简介多元线性回归多元线性回归模型实际中有很多问题是一个因变量与多个自变量成线性相关,我们可以用一个多元线性回归方程来表示。为了方便计算,我们将上式写成矩阵形式:Y=XW假设自变量维度为

多元线性回归

多元线性回归模型

实际中有很多问题是一个因变量与多个自变量成线性相关,我们可以用一个多元线性回归方程来表示。

为了方便计算,我们将上式写成矩阵形式:

Y = XW

假设自变量维度为NW为自变量的系数,下标0 - NX为自变量向量或矩阵,X维度为N,为了能和W0对应,X需要在第一行插入一个全是1的列。Y为因变量那么问题就转变成,已知样本X矩阵以及对应的因变量Y的值,求出满足方程的W,一般不存在一个W是整个样本都能满足方程,毕竟现实中的样本有很多噪声。最一般的求解W的方式是最小二乘法。

最小二乘法

我们希望求出的W是最接近线性方程的解的,最接近我们定义为残差平方和最小,残差的公式和残差平方和的公式如下:

上面的公式用最小残差平方和的方式导出的,还有一种思路用最大似然的方式也能推导出和这个一样的公式,首先对模型进行一些假设:

误差等方差不相干假设,即每个样本的误差期望为0,每个样本的误差方差都为相同值假设为σ误差密度函数为正态分布 e ~ N(0, σ^2)

简单推导如下:

由此利用最大似然原理导出了和最小二乘一样的公式。

最小二乘法求解

二次函数是个凸函数,极值点就是最小点。只需要求导数=0解出W即可。

模拟数据

我们这里用R语言模拟实践一下,由于我们使用的矩阵运算,这个公式一元和多元都是兼容的,我们为了可视化方便一点,我们就用R语言自带的women数据做一元线性回归,和多元线性回归的方式基本一样。women数据如下

> women height weight1 58 1152 59 1173 60 1204 61 1235 62 1266 63 1297 64 1328 65 1359 66 13910 67 14211 68 14612 69 15013 70 15414 71 15915 72 164

体重和身高具有线性关系,我们做一个散点图可以看出来:

我们用最小二乘推导出来的公式计算w如下

X <- cbind(rep(1, nrow(women)), women$height)X.T <- t(X)w <- solve(X.T %*% X) %*% X.T %*% y> w [,1][1,] -87.51667[2,] 3.45000> lm.result <- lm(women$weight~women$height)> lm.resultCall:lm(formula = women$weight ~ women$height)Coefficients: (Intercept) women$height -87.52 3.45

上面的R代码w使我们利用公式计算出来的,下边是R语言集成的线性回归函数拟合出来的,可以看出我们的计算结果是正确的,lm的只是小数点取了两位而已,将回归出来的函数画到图中看下回归的效果。画图对应的R代码如下,用R的感觉.....太飘逸了。

> png(file="chart2.png")> plot(women$height, women$weight)> lines(women$height, X %*% w)> dev.off()

梯度下降法

除了用正规方程方式求解W,也可以用最常见的梯度下降法求得W,因为最小二乘是个凸函数,所以这里找到的极小点就是最小点。下面这段代码用R写还是非常容易的,但是刚开始step步长参数调的太大了,导致一直不收敛,我还以为是程序错误,后来怎么看也没写错,就把参数调了个很小值,结果就收敛了。step的这个取值其实应该是变化的,先大后下比较科学,我这个调的很小,需要接近500万次才能收敛。

初始化W 为全0向量,也可以随机一个向量设置最大迭代次数,本例为了收敛设置了一个很大的数设置步长step,小了收敛很慢,大了不收敛.......求损失函数的梯度W(k+1) 为 W(k) + 损失函数负梯度 * 步长step循环,直到梯度接近0

X <- cbind(rep(1, nrow(women)), women$height)Y <- women$weightmaxIterNum <- 5000000;step <- 0.00003;W <- rep(0, ncol(X))for (i in 1:maxIterNum){ grad <- t(X) %*% (X %*% W - Y); if (sqrt(as.numeric(t(grad) %*% grad)) < 1e-3){ print(sprintf("iter times=%d", i)); break; } W <- W - grad * step;}print(W);

输出

[1] "iter times=4376771"

print(W);[,1][1,] -87.501509[2,] 3.449768

更多精彩文章 http://h2cloud.org/

, 1, 0, 9);

文章评论

Top