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)

Answer 2

I'd use calculus to calculate the expression for the maximum point. Differentiation will eliminate some of the constants in the equation, so the calculation is easier if you know what that max value needs to be.

If I recall correctly, a simple functions of 1 variable have a maximum at f'(x) = 0 and f''(x) < 0. Check me on this.

So if your function is f(x):

f(x) = a0 + a1*x + a2*x*x
f'(x) = a1 + 2*a2*x
f''(x) = 2*a2

Set the second equation equal to zero to get the stationary point, then put that value of x into the third one to find out if it's a max or min.