How to fit a quadratic model knowing the maximum in R?

Question

I have a data frame x and y, and I know the maximum of y. I want to fit this data to a quadratic model. How can I do it in R knowing the maximum? If I didn't know the maximum, I would fit it with lm(y~x + I(x^2)). Can anyone has an idea about this? Thanks in advance!

Answer 1

You have to minimize the sum of squares subject to the constraint; lm doesn't allow for constraints like this, so you have to use a generic optimization function, such as optim. Here's one way it could be done.

Make up some data. Here I'll say the known maximum is 50.

set.seed(5)
d <- data.frame(x=seq(-5, 5, len=51))
d$y <- 50 - 0.3*d$x^2 + rnorm(nrow(d))
M <- 50

Make a function to get the quadratic curve for points at x with given quadratic and linear coefficients and given maximum M. The calculus is straightforward; see duffymo's answer for details.

qM <- function(a, b, x, M) {
  c <- M - (3*b^2)/(4*a)
  a*x^2 + b*x + c
}

Make a function that get the sum of squares between a quadratic curve with given quadratic and linear coefficients and the data in d.

ff <- function(ab, d, M) {
  p <- qM(ab[1], ab[2], d$x, M)
  y <- d$y
  sum((p-y)^2)
}

Get the ordinary lm fit to use as starting values.

m0 <- lm(y ~ I(x^2) + x, data=d)
start <- coef(m0)[2:3]

Optimize the coefficients in the ff function.

o <- optim(start, ff, d=d, M=M)
o$par

Make a plot showing how the fit has a max at 50; the original lm fit doesn't.

plot(d)
xs <- seq(-5, 5, len=101)
lines(xs, predict(m0, newdata=data.frame(x=xs)), col="gray")
lines(xs, qM(o$par[1], o$par[2], xs, M))
abline(h=50, lty=3)