how to compute maximum of fitted polynomial?

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

how to compute maximum of fitted polynomial?

Joseph Clark
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    
______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: how to compute maximum of fitted polynomial?

Rui Barradas
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  
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: how to compute maximum of fitted polynomial?

Bert Gunter
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 = "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
>> ______________________________________________
>> [hidden email] 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.
>>
>
> ______________________________________________
> [hidden email] 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: how to compute maximum of fitted polynomial?

Hans W Borchers
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 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: how to compute maximum of fitted polynomial?

David Winsemius

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 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: how to compute maximum of fitted polynomial?

Hans W Borchers
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/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: how to compute maximum of fitted polynomial?

Joseph Clark
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/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.