前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >多元线性回归公式推导及R语言实现

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

作者头像
知然
发布2018-09-28 15:05:00
1.4K0
发布2018-09-28 15:05:00
举报
文章被收录于专栏:H2Cloud

多元线性回归

多元线性回归模型

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

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

Y = XW

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

最小二乘法

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

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

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

简单推导如下:

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

最小二乘法求解

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

模拟数据

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

代码语言:javascript
复制
> women
   height weight
1      58    115
2      59    117
3      60    120
4      61    123
5      62    126
6      63    129
7      64    132
8      65    135
9      66    139
10     67    142
11     68    146
12     69    150
13     70    154
14     71    159
15     72    164

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

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

代码语言:javascript
复制
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.result

Call:
lm(formula = women$weight ~ women$height)

Coefficients:
 (Intercept)  women$height  
      -87.52          3.45

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

画图对应的R代码如下,用R的感觉.....太飘逸了。

代码语言:javascript
复制
> 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
代码语言:javascript
复制
X <- cbind(rep(1, nrow(women)), women$height)
Y <- women$weight
maxIterNum <- 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

本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
原始发表:2018-08-28 ,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 多元线性回归
    • 多元线性回归模型
      • 最小二乘法
        • 最小二乘法求解
          • 模拟数据
            • 梯度下降法
              • 输出
          领券
          问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档