Fitting a quadratic through 5 points, goal is to find the maximum

Question

I have some physical experiments done at various locations. The locations produces a set of observations y for one value of x, the independent variable. In the end across a set of locations I have values in the following form

[y11, y21, y31, y41, y51, y61...] for one value of x, say x1

Then I repeat the experiment and get a new set of values

[y21, y22, y32, y42, y52, y62...] for a different value of x say x2

Ans so on.

In the end I have readings for y for 5 distinct values of x, [x1, x2, x3, x4, x5]

I wish to fit a quadratic to this data, and my main goal is to find the value of x, for which y is maximum. One way to do this is to define an average y for each x1 and fit a quadratic with the averaged out values. I know that fitting a quadratic to just 5 points is not a good idea. I am open to other ideas, that can help me solve the problem without directly fitting a functional relationship.

Some non-parametric idea for instance. One idea I have is to do some distribution analysis on values of y for x1, vs values of y for x2 and so on. This enables me to take all the y values for one x1, without averaging them. I would be open to other ideas and suggestions.

Answer 1

You are right that inferring a parabola from five points is too little data. But you have more than five points, namely all your measurements! Don't do any averaging (though this would already help), just fit the parabola to all your data. Model fitting doesn't mind multiple x values.

Let's do this in R. Some dummy data over five different x values:

set.seed(1)
xx <- rep(1:5,each=10)
yy <- -xx^2+6*xx-5+rnorm(length(xx),0,1)

Now we can fit the model. Note the I() to protect your square term, and note that this really presupposes homoskedastic errors:

model <- lm(yy~xx+I(xx^2))

Now we have a quadratic relationship. Some elementary calculus gives us the (estimated) x coordinate of maximum of the fitted parabola:

xx.max <- -coef(model)[2]/(2*coef(model)[3])

It's always a good idea to get an idea about the variability of our results. It may be possibly to derive a confidence interval for xx.max analytically, but the bootstrap is always easier, and given enough data (50 points should be enough), it should be valid:

require(boot)
foo <- boot(data=data.frame(xx=xx,yy=yy),statistic=function(data,indices){
        model <- lm(yy~xx+I(xx^2),data[indices,])
        -coef(model)[2]/(2*coef(model)[3])},
    strata=xx,
    R=1000)

Note that I am doing a stratified bootstrap, i.e., I am sampling with replacement within each x value, which makes sense here, since the data really are stratified.

So we can plot our points, the fitted parabola and the x coordinate of the estimated maximum together with the bootstrapped confidence interval:

plot(xx,yy,pch=19)
rect(quantile(foo$t,0.025),-2,quantile(foo$t,0.975),6,border=NA,col="lightgrey")
points(xx,yy,pch=19)
xx.plot <- seq(1,5,by=.01)
lines(xx.plot,predict(model,newdata=data.frame(xx=xx.plot)))
abline(v=xx.max)