

My script fits a thirdorder 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
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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 = "LBFGSB", 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 04062013 21:32, Joseph Clark escreveu:
> My script fits a thirdorder 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
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
>
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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 < [hidden email]> 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 = "LBFGSB", 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 04062013 21:32, Joseph Clark escreveu:
>>
>> My script fits a thirdorder 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
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/rhelp>> PLEASE do read the posting guide
>> http://www.Rproject.org/postingguide.html>> and provide commented, minimal, selfcontained, reproducible code.
>>
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.

Bert Gunter
Genentech Nonclinical Biostatistics
Internal Contact Info:
Phone: 4677374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdbfunctionalgroups/pdbbiostatistics/pdbncbhome.htm______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Bert Gunter <gunter.berton <at> 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 04062013 21:32, Joseph Clark escreveu:
> >>
> >> My script fits a thirdorder 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
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


On Jun 4, 2013, at 10:15 PM, Hans W Borchers wrote:
> Bert Gunter <gunter.berton <at> 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 04062013 21:32, Joseph Clark escreveu:
>>>>
>>>> My script fits a thirdorder 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
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


David Winsemius <dwinsemius <at> 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
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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.
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.

