[R] solving x in a polynomial function

[Mike Rennie]

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