我想计算参数的方差-协方差矩阵。参数通过非线性最小二乘拟合获得。
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)
现在尝试用粗麻布来做这件事,我无法让它工作以处理下面的错误消息
粗麻布
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()
工作。
我在这里的数学能力不是很强,因此采用了反复试验的方法。
最佳答案
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 的非线性回归分析及其应用程序。约翰福克斯也有一个短online appendix您可能会发现有帮助。
关于r - 在 R 中计算向量值 Hessian,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/18796529/