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