Skip to content
Prev 318553 / 398503 Next

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: