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

Tag: r , max , regression , fit , quadratic Author: wbtwang Date: 2011-10-11

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!

Do you really mean you know the maximum of y or you know the x coordinate of the maximum of y?
Would like to know the answer to the previous comment. If the answer is "maximum of y" then @Aaron's elegant answer works. If the answer is "x coordinate" then lm(y~I((x-xcoordmax)^2)) should work.
Kinda makes you think twice about answering questions from userXXXXXX people with near-zero rep. Aaron deserves mega-points for that answer.
Oh, I don't know. Other people seem to have found it useful, even if the OP's disappeared. Plus it's like a five-minute answer if you've done anything like it before. But thanks for the kudos, @Spacedman.
Sorry guys,I was out of town that is why I was not following your answers and comments. Thank you very much for your cooperation.

Best Answer

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)

image comparing lm fit with my fit

Other Answer1

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.