2

我试图从我自己的数据中获取 Hessian 矩阵，我有两个结果 -

  • 使用库中的代码 Hessian(numDeriv)
  • 使用库中的代码 numericHessian (maxLik)

相对于 numericHessian 的结果，Hessian 的结果非常小。

在这种情况下，我应该相信哪些结果？

具体来说，我使用的数据范围从 350000 到 1100000，它们是 9X2 矩阵，共有 18 个数据值。

我使用了一种标准偏差公式，“numericHessian”的结果范围从 230 到 466，使用 2X2 矩阵，而“Hessian”的结果范围从 -3.42e-18 到 1.34e-17，远远小于上一个。

您认为哪种标准差的计算正确？

代码如下：

 data=read.table("C:/file.txt", header=T);
 data <- as.matrix(data);
 library(plyr)
 library(MASS)
 w1 = tail(data/(rowSums(data)),1)
 w2 = t(w1)
 f <- function(x){
 w1 = tail(x/(rowSums(x)),1)
 w2 = t(w1)
 r = ((w1%*%cov(cbind(x))%*%w2)^(1/2))
 return(r)
 }
 library(maxLik); 
 numericHessian(f, t0=rbind(data[1,1], data[1,2]))
 library(numDeriv);
 hessian(f, rbind(data[1,1], data[1,2]), method="Richardson")

file.txt 如下：

  1                  2

 137                201

 122                342

 142                111

 171                126

 134                123

 823                876

 634                135

 541                214

 423                142

“numericHessian”的结果是：

          [,1]        [,2]
 [1,] 0.007105427 0.007105427
 [2,] 0.007105427 0.000000000

那么，“Hessian”的结果是：

          [,1]          [,2]
 [1,] -3.217880e-15 -1.957243e-16
 [2,] -1.957243e-16  1.334057e-16

非常感谢您提前。

4

1 回答 1

4

你没有给出一个可重复的例子，但我还是会尝试。

library(bbmle)
 x <- 0:10
 y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
 d <- data.frame(x,y)
 LL <- function(ymax=15, xhalf=6)
     -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE))
 fit <- mle2(LL)
 cc <- coef(fit)

以下是 MLE 中负对数似然函数的 Hessians（二阶导数矩阵）的有限差分估计：对这些矩阵求逆可以估计参数的方差-协方差矩阵。

 library(numDeriv)
 hessian(LL,cc)
 ##               [,1]          [,2]
 ## [1,]  1.296717e-01 -1.185789e-15
 ## [2,] -1.185789e-15  4.922087e+00

 library(maxLik)
 numericHessian(LL, t0=cc)
 ##           [,1]     [,2]
 ## [1,] 0.1278977 0.000000
 ## [2,] 0.0000000 4.916956

所以对于这个相对琐碎的例子，numDeriv::hessian并maxLik::numericHessian给出非常相似的结果。所以一定有一些你没有向我们展示的东西，或者你的问题的数字有什么特别之处。为了进一步进行，我们需要一个可重现的示例，请...

dat <- matrix(c(137,201,122,342,142,111,
                171,126,134,123,823,876,
                634,135,541,214,423,142),
        byrow=TRUE,ncol=2)
f <- function(x){
    w1 <- tail(x/(rowSums(x)),1)
    sqrt(w1%*%cov(cbind(x))%*%t(w1))
}
p <- t(dat[1,1:2,drop=FALSE])
f(p)  ## 45.25483
numDeriv::hessian(f,p)
##              [,1]          [,2]
## [1,] -3.217880e-15 -1.957243e-16
## [2,] -1.957243e-16  1.334057e-16
maxLik::numericHessian(f,t0=p)
##            [,1]        [,2]
## [1,] 0.007105427 0.007105427
## [2,] 0.007105427 0.000000000

好的，这些显然不同意。我不知道为什么，但在这种特殊情况下，我们可以分析你在做什么，看看哪个是正确的：

  • 由于您的输入矩阵是单列，x/rowSums(x)是一个向量，所以最后一个元素 ( w1 <- tail(...,1)) 只是 1。
  • 所以你的表达减少到sqrt(cov(cbind(x)))。同样，因为x是单列矩阵，cov()只是方差，sqrt(cov(.))只是标准差，或向量的范数。
  • 方差是任何元素与平均值的偏差的二次函数，因此标准偏差或多或少与平均值的偏差呈线性关系（零除外），因此我们预计二阶导数为零。所以看起来numDeriv::hessian给出了正确的答案

我们也可以通过增加epsfor来确认这一点numericHessian：

maxLik::numericHessian(f,t0=p,eps=1e-3)
##      [,1]          [,2]
## [1,]    0  0.000000e+00
## [2,]    0 -7.105427e-09

底线是使用更准确（但更慢）的方法，但如果你小心的话numDeriv，你可以得到合理的答案。numericHessian

于 2014-04-05T14:16:36.410 回答