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 <ehl...@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 <smartpink...@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 <ruipbarra...@sapo.pt> >>> To: Mike Rennie <mikerenni...@gmail.com> >>> Cc: r-help Mailing List <r-help@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@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. >>> >> >> > > -- Michael Rennie, Research Scientist Fisheries and Oceans Canada, Freshwater Institute Winnipeg, Manitoba, CANADA ______________________________________________ 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.