How to estimate the uncertainty in the zeros of a fitted function?

Question

I have fitted points with a polynomial. I now have the coefficients and the covariance matrix.

For a given y (in this case y=0; that is, x is a root of the polynomial) what is the uncertainty of that x where y=f(x)?

Answer 1

The principal objective of this reply is to point out how perilous this enterprise can be. Along the way I'll be able to suggest some approaches as well as provide some ideas for a different analysis. Whether any of this works will depend on the details of your circumstances.

The key points to watch for are

1. You need to get the model right. In particular, polynomial regression likely is going to do a poor job. Use splines instead.

2. It will be difficult to quantify the uncertainty in zeros located near stationary (near-level) points of the function.

3. Simulation (equivalently, a parametric bootstrap) can reveal much with relatively little effort.

---

The issues are best explained with an illustration. I constructed a quartic polynomial function on the interval $[-1,1]$ that, when raised or lowered a little, can have anywhere from zero through four real zeros. Here is a graph of it in black with red dots showing its zeros:

The open circles form a dataset of 37 points sampled at even intervals across the domain. Their values include iid Gaussian noise with a standard deviation of $0.1.$

The main problem is that small chance differences in that noise can cause any fitted curve to miss two (or occasionally all four) of the zeros and, on relatively rare occasions, may cause two or three of them to merge. Thus, not only are the locations of the zeros uncertain, even their number may be uncertain.

To illustrate this, I generated 400 such datasets, fitted a quartic polynomial to each dataset, found its zeros, and plotted their locations as vertical lines:

The most positive zero near $0.7$ is consistently estimated, but the others are all over the place. I ran a cluster analysis of their locations, resulting in finding five apparent clusters. The colors distinguish them. Posted above the cluster centers are percentages: these are the proportions of the 400 datasets in which a zero was found in each cluster: they estimate the chance that you will even detect a zero within each cluster.

If this weren't amusing enough, notice the reference to "degree 4 fits." The zeros were found by fitting a degree-4 polynomial to each dataset and then numerically finding all the zeros of that polynomial, of which there can be no more than four. In real life we usually don't know the correct degree. What happens when we specify a degree that isn't the same as the underlying function (or cannot approximate it well)? Here's what happens to the same 400 synthetic datasets when using degree-3 fits:

This procedure consistently gets things totally wrong: it always finds exactly two zeros and they're almost always in the wrong places. The moral is that you must use a fitting procedure that is capable of reproducing the true underlying function. Polynomial regression usually doesn't do that unless you're lucky. Use a spline or some similarly flexible method instead.

Compared to the first diagram (the degree-4 fit), this procedure appears to do a better job at positioning the zeros and it gets the right number of them more frequently.

Finally, as these plots suggest, you can use the spread of zeros within each cluster to summarize part of the uncertainty in the locations of the zeros. As you can see from this last plot, though, that's only part of the picture: how do you interpret five clusters of zeros when it's pretty clear the underlying function likely has only four zeros? Perhaps you should be content with a graphical illustration like this rather than with summary statistics. How you choose to summarize the uncertainty in the number of zeros will depend on your application and the interpretation of those zeros.

If you're lucky all your zeros will be "strongly transverse" ones like the high zero near $0.7$ in this pictures. In such cases, the standard error of the location of the zero will be proportional to the residual standard error and inversely proportional to the slope of the fitted curve at that zero. (There will not be a universal constant of proportionality, though.)

---

As requested, here is the R code used to generate the figures. It includes functions for finding all zeros and a section that performs a (quick and dirty) cluster analysis of a set of zeros.

#
# Find a zero of a function described by parallel arrays (x,y) where `x` is
# sorted in increasing order.
#
zero <- function(x, y, ...) {
  if (prod(range(y)) > 0) return(c())    # No zeros exist

  j <- min(which(y[-1] * y[1] < 0)) + 1  # Search for a change of sign from y[1]
  i <- 1
  while (y[i] * y[j] <= 0) i <- i+1      # Find the point just before the change
  i <- max(1,i-1)
  j <- min(length(y),j)
  if (i==j) return(x[i])

  f <- splinefun(x[c(i,j)], y[c(i,j)])   # Interpolate to find the zero
  uniroot(f, x[c(i,j)], ...)$root
}
#
# Repeatedly call `zero` to find all zeros.
#
zeros <- function(x, y, depth=0, tol=1e-4, ...) {
  if (depth >= 10) return(c())           # Avoids stack overflow
  tol.this <- tol * diff(range(x))
  x.0 <- zero(x, y, ...)

  # Recursively find zeros to the left and right of `x.0`:
  x.l <- x.u <- c()
  if (!is.null(x.0)) {
    l <- x <= x.0 - tol.this
    u <- x >= x.0 + tol.this
    if (sum(l) > 1) x.l <- zeros(x[l], y[l], depth+1, tol, ...)
    if (sum(u) > 1) x.u <- zeros(x[u], y[u], depth+1, tol, ...)
  }
  c(x.l, x.0, x.u) # This keeps the zeros in ascending order
}
#------------------------------------------------------------------------------#
library(splines)
set.seed(17)

x <- seq(-1, 1, length.out=37)                   # Fixed regressors
beta <- c(-1/8 + 0.02 + 1/16, 1/9.8, 1, 0, -2)   # Polynomial coefficients
y.0 <- outer(x, 1:length(beta)-1, `^`) %*% beta  # True values
sigma <- 0.1                                     # Gaussian error SD
degree <- 4                                      # Degree (or DF) to fit
method <- c("Polynomial", "Spline")[2]           # Fitting method
#
# Pretending `beta` is an estimate from data, perform a parametric bootstrap
# to explore the distributions of zeros.
#
N <- 4e2                                                 # Number of replications
Y <- data.frame(x = seq(min(x), max(x), length.out=201)) # Predict values here

Z <- replicate(N, {
  X <- data.frame(x = x, y = y.0 + rnorm(length(y.0), 0, sigma))
  if (method=="Polynomial") {
    fit <- lm(y ~ poly(x, degree=degree), X)
  } else {
    fit <- lm(y ~ bs(x, df=degree), X)
  }
  zeros(Y$x, predict(fit, newdata=Y))
})
#
# Usually `Z` will be a list, but in case all its elements are the same length
# `replicate` converts it into a matrix.
#
if("list" %in% class(Z)) z <- unlist(Z) else z <- c(Z)
#
# Perform a cluster analysis.  For illustrative purposes this is done
# automatically; in practice it might be better to do it in a supervised,
# exploratory mode in order to learn more about the patterns of zeros.
#
h <- hclust(dist(z))
k <- 0                # Number of clusters of zeros
while (k < 10) {      # Search for a reasonable number of clusters
  k <- k+1
  g <- cutree(h, k=k)
  omega <- tabulate(g, max(g)) / N
  if (max(omega) <= 1) break
}
x.0 <- by(z, g, mean) # Estimate cluster centers
#
# Plot the results.
#
cols <- terrain.colors(length(omega)+2, alpha=1/8)[1:length(omega)]
X <- data.frame(x = x, y = y.0 + rnorm(length(y.0), 0, sigma)) # Example dataset
main <- if(method=="Polynomial") {
  paste("A polynomial function fit with a degree", degree, "polynomial")
} else {
  paste("A polynomial function fit by cubic spline with", degree, "d.f.")
}

plot(x, y.0, type="l", lwd=2, ylab="y", cex.main=1, main=main) # True graph
abline(h=0, col="Red", lwd=2)                                  # y=0
abline(v=z, col=cols[g])                                       # Boootstrap zeros
mtext(sprintf("%.0f%%", 100*omega), at=x.0, cex=0.9)           # Cluster proportions
points(u, rep(0, length(u)), pch=21, bg="Red")
u <- zeros(x, y.0)                                             # True zeros
# with(X, points(x, y))                                        # Example data