[R] solving x in a polynomial function
William Dunlap
wdunlap at tibco.com
Fri Mar 1 23:11:23 CET 2013
You can avoid doing the algebra by using uniroot() to numerically find where
the predicted value hits your desired value. E.g.,
> fit <- lm(a ~ poly(b, 4))
> invertFit <- function(a){
+ stopifnot(length(a)==1)
+ uniroot(function(b)predict(fit, newdata=list(b=b))-a, interval=range(b))$root
+ }
> invertFit(5)
[1] 3.506213
See that it works with a plot:
> plot(b, a)
> btmp <- seq(min(b), max(b), len=129)
> lines(btmp, predict(fit, newdata=list(b=btmp)))
> abline(h=5, v=invertFit(5))
> abline(h=7, v=invertFit(7))
uniroot will not tell you if there is a problem due to the fit being nonmontonic, so
check the plot. (E.g., try lm(a ~ poly(b, 8)).)
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of Mike Rennie
> Sent: Friday, March 01, 2013 1:48 PM
> To: Peter Ehlers
> Cc: R help
> Subject: Re: [R] solving x in a polynomial function
>
> Hi Peter,
>
> With the edit you suggested, now I'm just getting back the value of a
> that I put in, not the expected value of b...
>
> > cbind(a,b)
> a b
> [1,] 1 1.0
> [2,] 2 2.0
> [3,] 3 2.5
> [4,] 4 3.0
> [5,] 5 3.5
> [6,] 6 4.0
> [7,] 7 6.0
> [8,] 8 7.0
> [9,] 9 7.5
> [10,] 10 8.0
>
> #a of 5 should return b of 3.5
>
> > realroots <- function(model, b){
> + is.zero <- function(x, tol = .Machine$double.eps^0.5) abs(x) < tol
> + if(names(coef(model))[1] == "(Intercept)")
> +
> + r <- polyroot(c(coef(model)[1] - b, coef(model)[-1]))
> + else
> + r <- polyroot(c(-b, coef(model)))
> + Re(r[is.zero(Im(r))])
> + }
> >
> > r <- realroots(po.lm, 5)
> > predict(po.lm, newdata = data.frame(b = r))
> 1 2
> 5 5
>
>
> This function just returns what I feed it as written.
>
> Mike
>
> On 3/1/13, Peter Ehlers <ehlers at ucalgary.ca> wrote:
> > On 2013-03-01 13:06, Mike Rennie wrote:
> >> Hi guys,
> >>
> >> Perhaps on the right track? However, not sure if it's correct. I fixed
> >> the bug that A.K. indicated above (== vs =), but the values don't seem
> >> right. From the original data, an a of 3 should give b of 2.5.
> >>
> >>> realroots <- function(model, b){
> >> + is.zero <- function(x, tol = .Machine$double.eps^0.5) abs(x) <
> >> tol
> >> + if(names(model)[1] == "(Intercept)")
> >> + r <- polyroot(c(coef(model)[1] - b, coef(model)[-1]))
> >> + else
> >> + r <- polyroot(c(-b, coef(model)))
> >> + Re(r[is.zero(Im(r))])
> >> + }
> >>>
> >>> r <- realroots(po.lm, 3)
> >>> predict(po.lm, newdata = data.frame(b = r)) # confirm
> >> 1
> >> 1.69
> >>
> >> So I think there's a calculation error somehwere.
> >
> > You need to replace the following line
> >
> > if(names(model)[1] == "(Intercept)")
> >
> > with
> >
> > if(names(coef(model))[1] == "(Intercept)")
> >
> >
> > Peter Ehlers
> >
> >
> >>
> >>
> >> On 3/1/13, arun <smartpink111 at yahoo.com> wrote:
> >>> Hi Rui,
> >>>
> >>> Looks like a bug:
> >>> ###if(names(model)[1] = "(Intercept)")
> >>> Should this be:
> >>>
> >>> if(names(model)[1] == "(Intercept)")
> >>>
> >>> A.K.
> >>>
> >>>
> >>>
> >>> ----- Original Message -----
> >>> From: Rui Barradas <ruipbarradas at sapo.pt>
> >>> To: Mike Rennie <mikerennie.r at gmail.com>
> >>> Cc: r-help Mailing List <r-help at r-project.org>
> >>> Sent: Friday, March 1, 2013 3:18 PM
> >>> Subject: Re: [R] solving x in a polynomial function
> >>>
> >>> Hello,
> >>>
> >>> Try the following.
> >>>
> >>>
> >>> a <- 1:10
> >>> b <- c(1, 2, 2.5, 3, 3.5, 4, 6, 7, 7.5, 8)
> >>>
> >>> dat <- data.frame(a = a, b = b) # for lm(), it's better to use a df
> >>> po.lm <- lm(a~b+I(b^2)+I(b^3)+I(b^4), data = dat); summary(po.lm)
> >>>
> >>>
> >>> realroots <- function(model, b){
> >>> is.zero <- function(x, tol = .Machine$double.eps^0.5) abs(x) < tol
> >>> if(names(model)[1] = "(Intercept)")
> >>> r <- polyroot(c(coef(model)[1] - b, coef(model)[-1]))
> >>> else
> >>> r <- polyroot(c(-b, coef(model)))
> >>> Re(r[is.zero(Im(r))])
> >>> }
> >>>
> >>> r <- realroots(po.lm, 5.5)
> >>> predict(po.lm, newdata = data.frame(b = r)) # confirm
> >>>
> >>>
> >>>
> >>> Hope this helps,
> >>>
> >>> Rui Barradas
> >>>
> >>> Em 01-03-2013 18:47, Mike Rennie escreveu:
> >>>> Hi there,
> >>>>
> >>>> Does anyone know how I solve for x from a given y in a polynomial
> >>>> function? Here's some example code:
> >>>>
> >>>> ##example file
> >>>>
> >>>> a<-1:10
> >>>>
> >>>> b<-c(1,2,2.5,3,3.5,4,6,7,7.5,8)
> >>>>
> >>>> po.lm<-lm(a~b+I(b^2)+I(b^3)+I(b^4)); summary(po.lm)
> >>>>
> >>>> (please ignore that the model is severely overfit- that's not the
> >>>> point).
> >>>>
> >>>> Let's say I want to solve for the value b where a = 5.5.
> >>>>
> >>>> Any thoughts? I did come across the polynom package, but I don't think
> >>>> that does it- I suspect the answer is simpler than I am making it out
> >>>> to be. Any help would be welcome.
> >>>>
> >>>
> >>> ______________________________________________
> >>> R-help at 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 at 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.
> >>>
> >>
> >>
> >
> >
>
>
> --
> Michael Rennie, Research Scientist
> Fisheries and Oceans Canada, Freshwater Institute
> Winnipeg, Manitoba, CANADA
>
> ______________________________________________
> R-help at 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.
More information about the R-help
mailing list