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

时间：2021-08-14 18:07:58

多元线性回归

多元线性回归模型

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

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

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

Y = XW

  • 假设自变量维度为N
  • W为自变量的系数，下标0 - N
  • X为自变量向量或矩阵，X维度为N,为了能和W0对应，X需要在第一行插入一个全是1的列。
  • Y为因变量

    那么问题就转变成，已知样本X矩阵以及对应的因变量Y的值，求出满足方程的W，一般不存在一个W是整个样本都能满足方程，毕竟现实中的样本有很多噪声。最一般的求解W的方式是最小二乘法。

最小二乘法

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

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

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

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

简单推导如下:

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

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

最小二乘法求解

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

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

模拟数据

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

women数据如下

> 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

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

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

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

X <- cbind(rep(1, nrow(women)), women$height)

X.T <- t(X)

Y <- women$weight

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代码如下，用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

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

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

归一化

上面的批量梯度下降为什么收敛如此之慢呢？原因很简单，没有做归一化，做了归一化，收敛速度快了非常非常多！！！！

正确代码如下：

XScale = scale(women$height)

Ux = attr(XScale, "scaled:center")

Dx = attr(XScale, "scaled:scale")

YScale = scale(women$weight)

Uy = attr(YScale, "scaled:center")

Dy = attr(YScale, "scaled:scale")

X <- cbind(rep(1, nrow(women)), as.matrix(XScale))

Y <- as.matrix(YScale)

maxIterNum <- 5000;

step <- 0.001;

W <- rep(0, ncol(X))

for (i in 1:maxIterNum){

    grad <- t(X) %*% (X %*% W -  Y);

    if (sqrt(as.numeric(t(grad) %*% grad)) < 1e-6){

        print(sprintf('iter times=%d', i));

        break;

    }

    W <- W - grad * step;

}

print(W);

W0 = W[1]

Wn = W[2:length(W)]

Wn = Dy * Wn / Dx

W0 = Uy + Dy * W0 - Dy * Ux / Dx

W  = c(W0, Wn)

print(W);

输出

[1] "iter times=1168"

print(W);

-88.53154 3.45000

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

相关文章