r - 计算R中的向量值Hessian

Question

我想计算参数的方差-协方差矩阵。通过非线性最小二乘拟合获得参数。

library(minpack.lm)
library(numDeriv)

变数

t <- seq(0.1,20,0.3)
a <- 20
b <- 14
c <- 0.4
jitter <- rnorm(length(t),0,0.5)
Hobs <- a+b*exp(-c*t)+jitter

功能定义

Hhat <- function(parList, t) {parList$a + parList$b*exp(-parL
Hhatde <- function(par, t) {par[1] + par[2]*exp(-par[3]*t)}st$c*t)}
residFun <- function(par, t, observed) observed - Hhat(par,t)

初始条件

parStart = list(a = 20, b = 10 ,c = 0.5)

nls.lm

library(minpack.lm)
out1 <- nls.lm(par = parStart, fn = residFun, observed = Hobs,
              t = t, control = nls.lm.control(nprint=0))

我希望手动计算通过 vcov(out1)返回的内容
我尝试使用：但 sigma和 vcov(out1)似乎并不相同

J <- jacobian(Hhatde, c(19.9508523,14.6586555,0.4066367 ), method="Richardson", 
method.args=list(),t=t)
sigma <- solve((t(J)%*%J))
vcov(out1)

现在试图用粗麻布做它，我不能让它为下面的错误消息cf工作

粗麻布的

H <- hessian(Hhatde, x = c(19.9508523,14.6586555,0.4066367 ), method="complex", method.args=list(),t=t)

Error in hessian.default(Hhatde, x = c(19.9508523, 14.6586555, 0.4066367),  : 
  Richardson method for hessian assumes a scalar valued function.

如何使我的 hessian()工作。

我在数学上不是很专一，因此我尝试反复尝试。

Answer 1

vcov(out1)返回模型中参数的比例方差-协方差矩阵的估计值。梯度的叉积的倒数solve(crossprod(J))返回未缩放的方差-协方差矩阵的估计值。比例因子是误差的估计方差。因此，要使用模型中的梯度和残差来计算缩放的方差-协方差矩阵（带有一些舍入误差）：

df = length(Hobs) - length(out1$par)                # degrees freedom
se_var = sum(out1$fvec^2) / df                      # estimated error variance
var_cov = se_var * solve(crossprod(J))              # scaled variance-covariance 

print(var_cov)                                      
print(vcov(out1))

要了解非线性回归和非线性最小二乘，您可能希望查看Seber＆Wild的非线性回归或Bates＆Watts的非线性回归分析及其应用。约翰·福克斯（John Fox）还有一个简短的 online appendix，您可能会有所帮助。