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
