多元线性回归

多元线性回归模型

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

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

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数据如下

> 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如下

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的感觉.....太飘逸了。

> 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$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/

多元线性回归公式推导及R语言实现的更多相关文章

  1. logistic逻辑回归公式推导及R语言实现

    Logistic逻辑回归 Logistic逻辑回归模型 线性回归模型简单,对于一些线性可分的场景还是简单易用的.Logistic逻辑回归也可以看成线性回归的变种,虽然名字带回归二字但实际上他主要用来二 ...

  2. 机器学习-线性回归(基于R语言)

    基本概念 利用线性的方法,模拟因变量与一个或多个自变量之间的关系.自变量是模型输入值,因变量是模型基于自变量的输出值. 因变量是自变量线性叠加和的结果. 线性回归模型背后的逻辑——最小二乘法计算线性系 ...

  3. R语言与概率统计(二) 假设检验

    > ####################5.2 > X<-c(159, 280, 101, 212, 224, 379, 179, 264, + 222, 362, 168, 2 ...

  4. R语言解读多元线性回归模型

    转载:http://blog.fens.me/r-multi-linear-regression/ 前言 本文接上一篇R语言解读一元线性回归模型.在许多生活和工作的实际问题中,影响因变量的因素可能不止 ...

  5. R 语言中的多元线性回归

    示例 sessionInfo() # 查询版本及系统和库等信息 # 工作目录设置 getwd() path <- "E:/RSpace" setwd(path) rm(lis ...

  6. 【R】多元线性回归

    R中的线性回归函数比较简单,就是lm(),比较复杂的是对线性模型的诊断和调整.这里结合Statistical Learning和杜克大学的Data Analysis and Statistical I ...

  7. R语言解读一元线性回归模型

    转载自:http://blog.fens.me/r-linear-regression/ 前言 在我们的日常生活中,存在大量的具有相关性的事件,比如大气压和海拔高度,海拔越高大气压强越小:人的身高和体 ...

  8. 机器学习(一) 从一个R语言案例学线性回归

    写在前面的话 按照正常的顺序,本文应该先讲一些线性回归的基本概念,比如什么叫线性回归,线性回规的常用解法等.但既然本文名为<从一个R语言案例学会线性回归>,那就更重视如何使用R语言去解决线 ...

  9. 【数据分析】线性回归与逻辑回归(R语言实现)

    文章来源:公众号-智能化IT系统. 回归模型有多种,一般在数据分析中用的比较常用的有线性回归和逻辑回归.其描述的是一组因变量和自变量之间的关系,通过特定的方程来模拟.这么做的目的也是为了预测,但有时也 ...

随机推荐

  1. 在Excel中制作雷达图

    雷达图的作用 雷达图是专门用来进行多指标体系比较分析的专业图表.从雷达图中可以看出指标的实际值与参照值的偏离程度,从而为分析者提供有益的信息.雷达图一般用于成绩展示.效果对比量化.多维数据对比等等,只 ...

  2. java同一个类不同方法间的同步

    对象的方法中一旦加入synchronized修饰,则任何时刻只能有一个线程访问synchronized修饰的方法.假设有个数据对象拥有写方法与读方法,多线程环境中要想保证数据的安全,需对该对象的读写方 ...

  3. WinForm给控件加入hint文字

    本文代码主要是参考别人的,仅为个人记录,方面后续使用~ 效果图: 主要代码在一个Win32Utility类中,代码如下: public static class Win32Utility { [Dll ...

  4. ls 排序

    ls 排序 首先我们通过man 来看看 ls其中有几项; -S     sort by file size : 按大小降序 --sort=WORD sort by WORD instead of na ...

  5. Android学习笔记(广播机制)

    1.Android的广播机制介绍 收听收音机也是一种广播,在收音机中有很多个广播电台,每个广播电台播放的内容都不相同.接受广播时广播(发送方)并不在意我们(接收方)接收到广播时如何处理.好比我们收听交 ...

  6. POJ3264——Balanced Lineup(线段树)

    本文出自:http://blog.csdn.net/svitter 题意:在1~200,000个数中.取一段区间.然后在区间中找出最大的数和最小的数字.求这两个数字的差. 分析:按区间取值,非常明显使 ...

  7. R的数据图形

    R支持4种图形类型: base graphics, grid graphics, lattice graphics,  ggplot2. Base graphics是R的默认图形系统. 一.  基本图 ...

  8. npm缺少css-loader,/style-compiler,stylus-loader问题,npm没有权限无法全局更新问题【已解决】

    ERROR in ./node_modules/css-loader!./node_modules/vue-loader/lib/style-compiler?{"vue":tru ...

  9. 《Java大学教程》—第4章 方法的实现

    4.2~3 声明.实现.调用4.4 数据传递:实参.形参.返回值4.6 变量作用域:局部变量(区域内访问).全局变量4.7 重载:运算符重载.方法重载-->多态 1.答:P67方法(method ...

  10. 指数型生成函数(EGF)学习笔记

    之前,我们学习过如何使用生成函数来做一些组合问题(比如背包问题),但是它面对排列问题(有标号)的时候就束手无策了. 究其原因,是因为排列问题的递推式有一些系数(这个待会就知道了),所以我们可以修改一下 ...