dlmModReg 的最大似然估计

机器算法验证 r 回归 最大似然 dlm
2022-04-07 22:32:06

我正在研究 R 包 dlm。到目前为止，它似乎是非常强大和灵活的包，具有良好的编程接口和良好的文档。

我已经能够成功地使用 dlmMLE 和 dlmModARMA 来估计 AR(1) 过程的参数：

u <- arima.sim(list(ar = 0.3), 100)
fit <- dlmMLE(u, parm = c(0.5, sd(u)),
              build = function(x)
                dlmModARMA(ar = x[1], sigma2 = x[2]^2))
fit$par

现在我正在尝试使用类似的代码来估计简单线性回归模型的参数：

r <- rnorm(100)
u <- -1*r + 0.5*rnorm(100)
fit <- dlmMLE(u, parm = c(0, 1),
              build = function(x)
                dlmModReg(x[1]*r, FALSE, dV = x[2]^2))
fit$par

我希望 fit$par 接近 c(-1, 0.5)，但我不断得到类似的东西

[1] -0.0002118851  0.4884367070

系数 -1 估计不正确。然而，奇怪的是噪声的方差被正确返回。

我知道如果初始值不好，最大似然估计可能会失败，但我观察到 dlmLL 返回的似然函数在第一个坐标中非常平坦。

所以我想知道：可以使用 dlm 来估计这样的模型吗？我相信该模型是“非奇异的”，但是我不确定如何在 dlm 中计算似然函数。

任何提示都非常感谢。

3个回答

我认为您的设置不正确。尝试这个：

set.seed(1234)
r <- rnorm(100)
X <- r
u <- -1*X + 0.5*rnorm(100)
MyModel <- function(x)  dlmModReg(X, FALSE, dV = x[1]^2)
fit <- dlmMLE(u, parm = c(0.3), build = MyModel)
mod <- MyModel(fit$par)
dlmFilter(u,mod)$a

您从 fit$par 的唯一元素中恢复对观测方差的估计：

> fit
$par
[1] 0.4431803

$value
[1] -20.69313

$counts
function gradient 
      17       17 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

而您对系数的估计（在您的情况下应该在 -1 左右）可以作为 的最后一个元素获得dlmFilter(u,mod)$a，它给出了处理新观察时的状态值：

    > dlmFilter(u,mod)$m
  [1]  0.0000000 -1.1486921 -1.2123431 -1.1172783 -1.1231454 -1.1170222
  [7] -1.0974931 -1.1377114 -1.0378758 -1.0927136 -1.0955372 -1.0120210
 [13] -0.9874791 -1.0036429 -1.0765513 -1.0678725 -1.0795124 -1.1568597
 [19] -1.2044821 -1.2056687 -1.2102896 -1.2938958 -1.2922945 -1.2670604
 [25] -1.1789594 -1.1570172 -1.1601590 -1.1417200 -1.1585501 -1.1608675
 [31] -1.1616278 -1.1744861 -1.1717561 -1.1715025 -1.1568086 -1.1451311
 [37] -1.1520867 -1.1379211 -1.1270897 -1.1048035 -1.1015793 -1.1054597
 [43] -1.0621750 -1.0621218 -1.0696813 -1.0807651 -1.0816893 -1.0647963
 [49] -1.0643440 -1.0667282 -1.0626404 -1.0623697 -1.0586265 -1.0571205
 [55] -1.0569135 -1.0579224 -1.0607623 -1.0582257 -1.0495232 -1.0494288
 [61] -1.0539632 -1.0555427 -1.0553468 -1.0491239 -1.0488604 -1.0491036
 [67] -1.0510551 -1.0576294 -1.0611296 -1.0628612 -1.0626451 -1.0573650
 [73] -1.0629577 -1.0647724 -1.0658052 -1.0823839 -1.0753808 -1.0747229
 [79] -1.0747762 -1.0615243 -1.0630352 -1.0697431 -1.0666448 -1.0617227
 [85] -1.0585460 -1.0583981 -1.0563544 -1.0567715 -1.0544349 -1.0573228
 [91] -1.0588404 -1.0639155 -1.0625845 -1.0578004 -1.0571034 -1.0602645
 [97] -1.0604838 -1.0586019 -1.0580891 -1.0587096 -1.0577559

希望这可以帮助。

下面是实现我的解决方案和 Paramonov 的解决方案的代码（稍作修改：我已更改原始dlmFilter(u,mod)$a发布的答案 dlmFilter(u,mod)$m）。

library(dlm)
set.seed(1234)
reps      <- 100
MyEstimates <- YourEstimates <- matrix(0,reps,2)
for (i in (1:reps) ) {
X <- r <- rnorm(100)
u <- -1*r + 0.5*rnorm(100)
#
fit <- dlmMLE(u, parm = c(1, sd(u)),
              build = function(x)
                dlmModReg(r, FALSE, dV = x[2]^2,
                          m0 = x[1], C0 = matrix(0)))
YourEstimates[i,] <- fit$par
#
MyModel <- function(x)  dlmModReg(X, FALSE, dV = x[1]^2)
fit <- dlmMLE(u, parm = c(0.3), build = MyModel)
mod <- MyModel(fit$par)
MyEstimates[i,] <- c(dlmFilter(u,mod)$m[101],fit$par[1])
}

当我运行上面的代码时，这就是我得到的：

> summary(YourEstimates)
       V1                V2         
 Min.   :-9.5284   Min.   :-0.5747  
 1st Qu.:-1.4280   1st Qu.: 0.4710  
 Median :-0.9795   Median : 0.4937  
 Mean   :-0.9737   Mean   : 0.4369  
 3rd Qu.:-0.5636   3rd Qu.: 0.5215  
 Max.   : 4.5222   Max.   : 0.5980  
> summary(MyEstimates)
       V1                V2         
 Min.   :-1.1099   Min.   :-0.6010  
 1st Qu.:-1.0266   1st Qu.: 0.4736  
 Median :-0.9974   Median : 0.4961  
 Mean   :-0.9938   Mean   : 0.4469  
 3rd Qu.:-0.9635   3rd Qu.: 0.5158  
 Max.   :-0.8390   Max.   : 0.5776

虽然第一组估计值对第二个参数给出了类似的估计值，但它偶尔会给出第一个参数的值。我认为原因是将状态“绑定”到其初始值

C0=matrix(0)

导致数值不稳定，但我不确定。无论如何，您可能想查看该问题。

在阅读了 dlmFilter 的帮助之后，我可以想出以下代码：

r <- rnorm(100)
u <- -1*r + 0.5*rnorm(100)
fit <- dlmMLE(u, parm = c(1, sd(u)),
              build = function(x)
                dlmModReg(r, FALSE, dV = x[2]^2,
                          m0 = x[1], C0 = matrix(0)))
fit$par

[1] -1.1330088  0.4788357
其它你可能感兴趣的问题
上一篇R 的 cmdscale() 的大 O 缩放 下一篇“这就是她所说的”问题的朴素贝叶斯分类