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