Doh- I'm a moron. I get it now. The last line is the confirmation that
function "realroots" is working. Sorry- late in the day on a friday.
Thanks everyone for your help with this- Uber-useful, and much appreciated.
Mike
On 3/1/13, Mike Rennie <mikerenni...@gmail.com> wrote:
> 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
>
--
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.
