R: Confidence intervals on non-linear fit with a non-analytic model

Question

I need to fit x-y data with a model, which is non-analytic. I have a function f(x) that calculates the model for each x numerically, but there is no analytical equation. For the fit, I use optim in R. I minimise RMS between the model and the data. It works well and returns reasonable parameters.

I would like to find confidence intervals (or at least standard errors) on the best-fitting parameters. I found on internet that this can be done from the Hessian matrix, but only if maximising log-likelihood function. I don't know how to do this, all I have is x, y and f(x) from which I find RMS. Alas, I have no good way of estimating errors on y.

How can I find confidence intervals on my fit parameters?

Edit: perhaps an example in R might help explaining what I'm asking for. This example uses a simple analytic function to fit the data, in my real case the function is non-analytic, so I cannot use, e.g., 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]

The best fitting parameters are a = 0.50 and b = 0.20. I need to find 95% confidence intervals on these.

Answer 1

This is a job for the bootstrap:

(1) create a large number of synthetic datasets x*. These are created by sampling from x with replacement the same number of data as were in x. For example, if your data is (1,2,3,4,5,6) an x* might be (5,2,4,4,2,3) (note that values might appear multiple times, or not at all because we are sampling with replacement)

(2) For each x*, calculate f(x*). If there are other parameters which don't depend on the data, don't change them. (so f(x,a,b,c) becomes f(x*,a,b,c) as long as a,b,c don't depend on x. Call these quantities f*.

(3) You can estimate anything you want from these f*. If you want the standard deviation of f(x), take the standard deviation of f*. If you want the 95% confidence interval, take the range from the 2.5 to the 97.5 percentiles of f*. More formally, if you want to estimate g(f(x)) you estimate it as g(f(x*)).

I should say this is a very practically-oriented explanation of the bootstrap. I have glossed over many theoretical details, but the bootstrap is near-universally applicable (basically as long as the thing you are trying to estimate actually exists, you are usually okay).

To apply this to the example you have given in your code:

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))