Re: [R] how to compute maximum of fitted polynomial?
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] 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: predict(p2,xmax) 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. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to compute maximum of fitted polynomial?
David Winsemius comcast.net> writes: > [...] > > On Jun 4, 2013, at 10:15 PM, Hans W Borchers wrote: > > > In the case of polynomials, "elementary math ... methods" can > > actually be > > executed with R: library(polynomial) # -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] 1.42265 > These look like the functions present in the 'polynom' package > authored by Bill Venables [aut] (S original), Kurt Hornik [aut, cre] > (R port), Martin Maechler. I wasn't able to find a 'polynomial' > package on CRAN. The 'mpoly' package by David Kahle offers > multivariate symbolic operations as well. > Sorry, yes of course, it should be `library(polynom)`. Somehow I'm making this mistake again and again. And one has to be a bit careful about complex roots. Hans Werner __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to compute maximum of fitted polynomial?
On Jun 4, 2013, at 10:15 PM, Hans W Borchers wrote: Bert Gunter gene.com> writes: 1. This looks like a homework question. We should not do homework here. 2. optim() will only approximate the max. 3. optim() is not the right numerical tool for this anyway. optimize() is. 4. There is never a guarantee numerical methods will find the max. 5. This can (and should?) be done exactly using elementary math rather than numerical methods. Cheers, Bert In the case of polynomials, "elementary math ... methods" can actually be executed with R: library(polynomial) # -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] 1.42265 Obviously, the same procedure will work for polynomials p0 of higher orders. These look like the functions present in the 'polynom' package authored by Bill Venables [aut] (S original), Kurt Hornik [aut, cre] (R port), Martin Maechler. I wasn't able to find a 'polynomial' package on CRAN. The 'mpoly' package by David Kahle offers multivariate symbolic operations as well. -- David. Hans Werner Em 04-06-2013 21:32, Joseph Clark escreveu: My script fits a third-order polynomial to my data with something like this: model <- lm( y ~ poly(x, 3) ) What I'd like to do is find the theoretical maximum of the polynomial (i.e. the x at which "model" predicts the highest y). Specifically, I'd like to predict the maximum between 0 <= x <= 1. What's the best way to accomplish that in R? Bonus question: can R give me the derivative or 2nd derivative of the polynomial? I'd like to be able to compute these at that maximum point. Thanks in advance! -- David Winsemius, MD Alameda, CA, USA __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to compute maximum of fitted polynomial?
Bert Gunter gene.com> writes: > > 1. This looks like a homework question. We should not do homework here. > 2. optim() will only approximate the max. > 3. optim() is not the right numerical tool for this anyway. optimize() is. > 4. There is never a guarantee numerical methods will find the max. > 5. This can (and should?) be done exactly using elementary math rather > than numerical methods. > > Cheers, > Bert In the case of polynomials, "elementary math ... methods" can actually be executed with R: library(polynomial) # -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] 1.42265 Obviously, the same procedure will work for polynomials p0 of higher orders. Hans Werner > > Em 04-06-2013 21:32, Joseph Clark escreveu: > >> > >> My script fits a third-order polynomial to my data with something like > >> this: > >> > >> model <- lm( y ~ poly(x, 3) ) > >> > >> What I'd like to do is find the theoretical maximum of the polynomial > >> (i.e. the x at which "model" predicts the highest y). Specifically, I'd > >> like to predict the maximum between 0 <= x <= 1. > >> > >> What's the best way to accomplish that in R? > >> > >> Bonus question: can R give me the derivative or 2nd derivative of the > >> polynomial? I'd like to be able to compute these at that maximum point. > >> > >> Thanks in advance! > >> > >> > >> // joseph w. clark , phd , visiting research associate > >> \\ university of nebraska at omaha - college of IS&T __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to compute maximum of fitted polynomial?
1. This looks like a homework question. We should not do homework here. 2. optim() will only approximate the max. 3. optim() is not the right numerical tool for this anyway. optimize() is. 4. There is never a guarantee numerical methods will find the max. 5. This can (and should?) be done exactly using elementary math rather than numerical methods. Cheers, Bert On Tue, Jun 4, 2013 at 2:51 PM, Rui Barradas wrote: > Hello, > > As for the first question, you can use ?optim to compute the maximum of a > function. Note that by default optim minimizes, to maximize you must set the > parameter control$fnscale to a negative value. > > fit <- lm(y ~ poly(x, 3)) > > fn <- function(x, coefs) as.numeric(c(1, x, x^2, x^3) %*% coefs) > > sol <- optim(0, fn, gr = NULL, coef(fit), control = list(fnscale = -1), > method = "L-BFGS-B", lower = 0, upper = 1) > > > As for the second question, I believe you can do something like > > dfdx <- D( expression(a + b*x + c*x^2 + d*x^3), "x") > > a <- coef(fit)[1] > b <- coef(fit)[2] > c <- coef(fit)[3] > d <- coef(fit)[4] > x <- sol$par > eval(dfdx) > > > See the help page for ?D > > > Hope this helps, > > Rui Barradas > > Em 04-06-2013 21:32, Joseph Clark escreveu: >> >> My script fits a third-order polynomial to my data with something like >> this: >> >> model <- lm( y ~ poly(x, 3) ) >> >> What I'd like to do is find the theoretical maximum of the polynomial >> (i.e. the x at which "model" predicts the highest y). Specifically, I'd >> like to predict the maximum between 0 <= x <= 1. >> >> What's the best way to accomplish that in R? >> >> Bonus question: can R give me the derivative or 2nd derivative of the >> polynomial? I'd like to be able to compute these at that maximum point. >> >> Thanks in advance! >> >> >> // joseph w. clark , phd , visiting research associate >> \\ university of nebraska at omaha - college of IS&T >> __ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to compute maximum of fitted polynomial?
Hello, As for the first question, you can use ?optim to compute the maximum of a function. Note that by default optim minimizes, to maximize you must set the parameter control$fnscale to a negative value. fit <- lm(y ~ poly(x, 3)) fn <- function(x, coefs) as.numeric(c(1, x, x^2, x^3) %*% coefs) sol <- optim(0, fn, gr = NULL, coef(fit), control = list(fnscale = -1), method = "L-BFGS-B", lower = 0, upper = 1) As for the second question, I believe you can do something like dfdx <- D( expression(a + b*x + c*x^2 + d*x^3), "x") a <- coef(fit)[1] b <- coef(fit)[2] c <- coef(fit)[3] d <- coef(fit)[4] x <- sol$par eval(dfdx) See the help page for ?D Hope this helps, Rui Barradas Em 04-06-2013 21:32, Joseph Clark escreveu: My script fits a third-order polynomial to my data with something like this: model <- lm( y ~ poly(x, 3) ) What I'd like to do is find the theoretical maximum of the polynomial (i.e. the x at which "model" predicts the highest y). Specifically, I'd like to predict the maximum between 0 <= x <= 1. What's the best way to accomplish that in R? Bonus question: can R give me the derivative or 2nd derivative of the polynomial? I'd like to be able to compute these at that maximum point. Thanks in advance! // joseph w. clark , phd , visiting research associate \\ university of nebraska at omaha - college of IS&T __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to compute maximum of fitted polynomial?
My script fits a third-order polynomial to my data with something like this: model <- lm( y ~ poly(x, 3) ) What I'd like to do is find the theoretical maximum of the polynomial (i.e. the x at which "model" predicts the highest y). Specifically, I'd like to predict the maximum between 0 <= x <= 1. What's the best way to accomplish that in R? Bonus question: can R give me the derivative or 2nd derivative of the polynomial? I'd like to be able to compute these at that maximum point. Thanks in advance! // joseph w. clark , phd , visiting research associate \\ university of nebraska at omaha - college of IS&T __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
