R：使用非分析模型进行非线性拟合的置信区间

Question

我需要用非分析模型拟合 x-y 数据。我有一个函数f(x) 以数值方式计算每个x 的模型，但没有解析方程。为了拟合，我在 R 中使用optim。我最小化模型和数据之间的 RMS。它运行良好并返回合理的参数。

我想找到最佳拟合参数的置信区间（或至少是标准误）。我在互联网上发现这可以从 Hessian 矩阵中完成，但前提是最大化对数似然函数。我不知道该怎么做，我只有x、y 和f(x) 从中找到RMS。唉，我没有很好的方法来估计 y 上的错误。

如何找到拟合参数的置信区间？

编辑：也许 R 中的一个示例可能有助于解释我的要求。这个例子使用了一个简单的解析函数来拟合数据，在我的真实案例中这个函数是非解析的，所以我不能使用，例如，nls。

set.seed(666)

# generate data
x <- seq(100) / 100
y <- 0.5 * x + rnorm(100, sd = 0.03) + 0.2

# function to fit
f <- function(x, a, b) {
  a * x + b
}

# error function to minimise: RMS
errfun <- function(par, x, y) {
  a <- par[1]
  b <- par[2]
  err <- sqrt(sum((f(x, a, b) - y)^2))
}

# use optim to fit the model to the data
par <- c(1, 0)
res <- optim(par, errfun, gr=NULL, x, y)

# best-fitting parameters
best_a <- res$par[1]
best_b <- res$par[2]

最佳拟合参数是 a = 0.50 和 b = 0.20。我需要找到这些的 95% 置信区间。

Answer 1

这是引导程序的工作：

(1) 创建大量合成数据集x*。这些是通过从x 采样创建的，并替换与x 中相同数量的数据。例如，如果您的数据是 (1,2,3,4,5,6)，那么 x* 可能是 (5,2,4,4,2,3)（请注意，值可能会出现多次，或者根本不会出现，因为我们正在采样替换）

(2) 对于每个x*，计算f(x*)。如果还有其他不依赖于数据的参数，请不要更改它们。 （所以f(x,a,b,c) 变为f(x*,a,b,c)，只要a,b,c 不依赖于x。将这些数量称为f*。

(3) 你可以从这些f* 中估计任何你想要的东西。如果您想要f(x) 的标准差，请取f* 的标准差。如果您想要 95% 的置信区间，请取 f* 的 2.5 到 97.5 个百分位数的范围。更正式地说，如果你想估计g(f(x))，你估计它为g(f(x*))。

我应该说这是对引导程序非常实用的解释。我已经掩盖了许多理论细节，但引导程序几乎是普遍适用的（基本上只要您尝试估计的事物确实存在，通常就可以了）。

将此应用于您在代码中给出的示例：

x <- seq(100) / 100
y <- 0.5 * x + rnorm(100, sd = 0.03) + 0.2

# function to fit
f <- function(x, a, b) {
  a * x + b
}

# error function to minimise: RMS
errfun <- function(par, x, y) {
  a <- par[1]
  b <- par[2]
  err <- sqrt(sum((f(x, a, b) - y)^2))
}

# this is the part where we bootstrap
# use optim to fit the model to the data
best_a <- best_b <- numeric(10000)
for(i in 1:10000){
  j <- sample(100,replace=TRUE)
  x.boot <- x[j]; y.boot <- y[j]
par <- c(1, 0)
res <- optim(par, errfun, gr=NULL, x.boot, y.boot)

# best-fitting parameters
best_a[i] <- res$par[1]
best_b[i] <- res$par[2]
}
# now, we look at the *vector* best_a
# for example, if you want the standard deviation of a,
sd(best_a)
# or a 95% confidence interval for b,
quantile(best_b,c(0.025,0.975))