R语言-回归系数的极大似然估计

老师要求我们对回归方程中的回归系数进行极大似然估计，回归方程如下：

计算步骤如下：
步骤一：写出似然函数log(β)，其中的β为(β0,β1,β2)t(β_0, β_1, β_2)^t(β0,β1,β2)t
步骤二：解出log’(β)=0的根

由回归分析的知识，yiy_iyi ~ N(β0+β1x1i+β2x2i,1)N( \beta_0 +\beta_1x_{1i}+\beta_2x_{2i},1)N(β0+β1x1i+β2x2i,1)
因此可以写出似然函数：

取对数

其中C为常数

但是在写代码的时候不用将极大似然函数表示出来，只需要对log’(β)=0求根即可，这里用的是牛顿法求根。

首先来简单介绍一下牛顿算法：
这是一种通过不断迭代来求根的算法，先随便给一个初始的x值，再给出迭代原则，最终就能找到方程的根。例如图中是求g(x)=3x5−4x4+6x3+4x−4g(x)=3x^5-4x^4+6x^3+4x-4g(x)=3x5−4x4+6x3+4x−4的根，通过三次迭代就可以很接近了。

图中展现的是一维牛顿迭代，从几何上易知：

从而推出

注意在我们的优化算法中g(x)实际上是f(x)的导数，因此一维牛顿迭代的公式为

现在把它推广到高维牛顿迭代：

βn+1=βn−H−1Gβ_{n+1} = β_n - H^{-1} Gβn+1=βn−H−1G

其中的H为黑塞矩阵，即二阶导矩阵

G为梯度矩阵，即一阶导矩阵

用牛顿迭代时，我们应该设置一个迭代停止的条件，不然程序会在零点附近永远运行下去。这里列出了一些可行的条件：

我设置的条件是 βn+1−βn<0.0001\beta_{n+1}-\beta_n < 0.0001βn+1−βn<0.0001

整个过程用R语言来实现：

#生成数据集
#y是n*1的矩阵，b是[b0,b1,b2],x是n*3的矩阵
set.seed(2)
bbb0 = c(10,10,10)
y = 10 + 10*rnorm(10, 1, 1)
x1 = c(rep(1,10))
x2 = y*2 + rnorm(10,0,1)
x3 = y - 5 + rnorm(10,0,1)
x = t(rbind(x1, x2, x3))#梯度矩阵
G = matrix(nrow = length(y), ncol = length(bbb0))
gloglike = function(b,y,x){for (i in 1:length(y)) {ei = y[i] - x[i,]%*%b   #是一个数ei = as.vector(ei)gi = ei * x[i,]   #1*3的矩阵G[i,]=gi}return(colSums(G))
}#黑塞矩阵
H = matrix(0, nrow = length(bbb0), ncol = length(bbb0))
Hloglike = function(x){for (i in 1:length(y)) {xi = as.matrix(x[i,])h = -xi%*%t(xi)H = H + h#hi = as.vector(h)#H[i,]=hi}#out = matrix(colSums(H)) #return(out)return(H)
}#牛顿求根算法
newton = function(x,y,bbb0){HH = Hloglike(x)GG = gloglike(bbb0,y,x)delta = c(2,3,4)n = 1
while (delta[1]>0.0001 || delta[2]>0.0001 || delta[3]>0.0001) {bbb1 = bbb0 - solve(HH)%*%GGdelta = abs(bbb1 - bbb0)bbb0 = bbb1GG = gloglike(bbb0,y,x)n = n + 1
}result = list(bbb0,n)return(result)
}newton(x,y,bbb0)

最后结果为

x1 0.51256600
x2 0.46233262
x3 0.06147585

n = 3

x1,x2,x3 即为 β0,β1,β2β_0, β_1, β_2β0,β1,β2 的估计结果，迭代次数为3次

R语言-回归系数的极大似然估计相关推荐

R语言-点估计与极大似然估计
这里写自定义目录标题前言一.点估计二.极大似然估计提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档文章目录前言一.点估计二.极大似然估计前言提示:这里可以添加本文要 ...
机器学习笔记1.矩估计、极大似然估计。
1.矩估计 1.1矩估计思想: 矩估计是基于一种简单的"替换"思想,即用样本矩估计总体矩 1.2矩估计理论: 矩估计的理论依据就是基于大数定律的,大数定律语言化表述为:当总体的k阶 ...
R语言作业一：矩估计、极大似然估计、拟合、对数正态分布、泊松分布、负二项分布
一.矩估计.极大似然估计.拟合.对数正态分布 ##导入数据 setwd("C:/Users/chang/Documents/SRM-PA/R简介/上课练习数据集") healthe ...
贝叶斯网专题11：参数学习之极大似然估计
第一部分:贝叶斯网基础 1.1 信息论基础 1.2 贝叶斯网基本概念 1.3 变量独立性的图论分析第二部分:贝叶斯网推理 2.1 概率推理中的变量消元方法 2.2 团树传播算法 2.3 近似推理 2 ...
似然函数的意义与极大似然估计
什么是概率? 简单来说,概率是一个函数,定义域是样本空间,满足非负性,规范性,可列可加性. 严格的公理化定义如下: 概率可以做什么?统计又可以做什么? 什么是先验概率,后验概率,似然? 先验概率:根据 ...
一文看懂 “极大似然估计” 与 “最大后验估计” —— 极大似然估计篇
参考: 唐宇迪<人工智能数学基础>第8章 Richard O. Duda <模式分类>第三章白板机器学习 P2 - 频率派 vs 贝叶斯派频率学派还是贝叶斯学派?聊一聊机器 ...
统计学笔记1：截尾分布的矩估计与极大似然估计
截尾分布的矩估计与极大似然估计在参数估计中,我们通常喜欢用极大似然估计来估计一个参数,这样估计的参数通常具有良好的性质,但有时其并不那么容易求解.在参数估计中,矩估计的计算方法较为简易,但其结果的偏 ...
基于接收信号强度（RSS）的室内定位/无线传感器网络定位——极大似然估计ML/最小二乘估计WLS
基于接收信号强度(RSS)的室内定位/无线传感器网络定位--极大似然估计ML/最小二乘估计WLS 原创不易,路过的各位大佬请点个赞针对AOA,TOA,TDOA,RSS等室内定位.导航的探讨.技术支持 ...
基于到达时间（TOA）的室内定位（/无线传感器网络定位）——极大似然估计ML
基于到达时间(TOA)的室内定位(/无线传感器网络定位)--极大似然估计ML 原创不易,路过的各位大佬请点个赞针对AOA,TOA,TDOA,RSS等室内定位.导航的探讨.技术支持.==代码(有偿)= ...

R语言-回归系数的极大似然估计

R语言-回归系数的极大似然估计相关推荐

最新文章

热门文章