[R] how to compute maximum of fitted polynomial?
joeclark77 at hotmail.com
Wed Jun 5 16:11:30 CEST 2013
Thank you all! This approach, using the 'polynom' library, did the trick.
> library(polynom) # -6 + 11*x - 6*x^2 + x^3
> p0 <- polynomial(c(-6, 11, -6, 1)) # has zeros at 1, 2, and 3
> p1 <- deriv(p0); p2 <- deriv(p1) # first and second derivative
> xm <- solve(p1) # maxima and minima of p0
> xmax = xm[predict(p2, xm) < 0] # select the maxima
> xmax #  1.42265
With these tweaks:
In fitting the model to the data I had to use raw=TRUE:
model <- lm( y ~ poly(x, 3, raw=TRUE) )
Then I could generate p0 directly from my lm:
p0 <- polynomial( coef(model) )
And I answered my second question by using the obvious:
I don't know what I would have done if the optima weren't between x=0 and x=1, which was my constraint. In that case the "maximum" would have been one of the endpoints rather than a zero of p1. I suppose I could just have checked for it with some if/then code. Fortunately it didn't turn out to be an issue with my data.
And no, this wasn't a homework problem. I didn't do the math by hand because I needed to automate this process for several subsets of my data and several fitted models.
More information about the R-help