Any better (more efficient, built-in) ideas for computing
coef[1]+coef[2]*x+coef[3]*x^2+ ...
than
polynom <- function(coef,x) {
n <- length(coef)
sum(coef*apply(matrix(c(rep(x,n),seq(0,n-1)),ncol=2),1,function(z)z[1]^z[2]))
}
?
Ben
--
318 Carr Hall bolker at zoo.ufl.edu
Zoology Department, University of Florida http://www.zoo.ufl.edu/bolker
Box 118525 (ph) 352-392-5697
Gainesville, FL 32611-8525 (fax) 352-392-3704
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Ben Bolker <ben at zoo.ufl.edu> writes:> Any better (more efficient, built-in) ideas for computing > > coef[1]+coef[2]*x+coef[3]*x^2+ ... > > than > > polynom <- function(coef,x) { > n <- length(coef) > > sum(coef*apply(matrix(c(rep(x,n),seq(0,n-1)),ncol=2),1,function(z)z[1]^z[2])) > }Back when debugging polyroot() I used this one polyeval<-function(coef,x){v<-0;for (i in rev(coef)) v<-v*x+i; v} Can't vouch for the efficiency, but it vectorizes nicely in x. -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Ben Bolker wrote:> > Any better (more efficient, built-in) ideas for computing > > coef[1]+coef[2]*x+coef[3]*x^2+ ... > > than > > polynom <- function(coef,x) { > n <- length(coef) > > sum(coef*apply(matrix(c(rep(x,n),seq(0,n-1)),ncol=2),1,function(z)z[1]^z[2])) > } > > ?I think the function polynom2 <- function(x, coef) sum(coef * x^seq(0,length(coef)-1)) is a little bit simpler and quicker for high order polynoms. But if you want to be able to compute the function for vector-valued x, then you might need something like polynom3 <- function(x, coef) apply(coef * sapply(as.matrix(x), function(x) x^seq(0,length(coef)-1)), 2, sum) Then you can do: R> mycoef <- rnorm(10000) R> polynom(mycoef, 0.5) [1] -0.1422343 R> polynom2(0.5, mycoef) [1] -0.1422343 R> polynom2(0.75, mycoef) [1] -1.395595 R> polynom3(c(0.5, 0.75), mycoef) [1] -0.1422343 -1.3955947 Best, Z> Ben > -- > 318 Carr Hall bolker at zoo.ufl.edu > Zoology Department, University of Florida http://www.zoo.ufl.edu/bolker > Box 118525 (ph) 352-392-5697 > Gainesville, FL 32611-8525 (fax) 352-392-3704 > > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- > r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > Send "info", "help", or "[un]subscribe" > (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Wed, 9 Oct 2002, Ben Bolker wrote:> > Any better (more efficient, built-in) ideas for computing > > coef[1]+coef[2]*x+coef[3]*x^2+ ... > > than > > polynom <- function(coef,x) { > n <- length(coef) > > sum(coef*apply(matrix(c(rep(x,n),seq(0,n-1)),ncol=2),1,function(z)z[1]^z[2])) > } >Well if x is a scalar as in your example n<-length(coef)-1 sum(coef*x^(0:n)) seems simpler, or if x is a vector then n<-length(coef)-1 rowSums(coef*outer(0:n,x,"^")) or if it's a long enough vector that you don't want n copies of it rval<-coef[1] xn<-x n<-length(coef) for(i in 2:n){ rval<-rval+coef[i]*xn xn<-xn*x } Unlike lapply, which is actually faster than a for() loop, apply() is basically a clarity optimisation rather than a speed optimisation, and in this case I don't think it's clearer. -thomas -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Hello Ben, would the supplementary package 'polynom' of help to you? library(polynom) p <- polynomial(c(1,-1.11,0,0,0.18)) predict(p, 'your values for x') Rgds, Bernhard -----Original Message----- From: Thomas Lumley [mailto:tlumley at u.washington.edu] Sent: 10 October 2002 00:43 To: bolker at zoo.ufl.edu Cc: R help list Subject: Re: [R] polynomial On Wed, 9 Oct 2002, Ben Bolker wrote:> > Any better (more efficient, built-in) ideas for computing > > coef[1]+coef[2]*x+coef[3]*x^2+ ... > > than > > polynom <- function(coef,x) { > n <- length(coef) > >sum(coef*apply(matrix(c(rep(x,n),seq(0,n-1)),ncol=2),1,function(z)z[1]^z[2]) )> } >Well if x is a scalar as in your example n<-length(coef)-1 sum(coef*x^(0:n)) seems simpler, or if x is a vector then n<-length(coef)-1 rowSums(coef*outer(0:n,x,"^")) or if it's a long enough vector that you don't want n copies of it rval<-coef[1] xn<-x n<-length(coef) for(i in 2:n){ rval<-rval+coef[i]*xn xn<-xn*x } Unlike lapply, which is actually faster than a for() loop, apply() is basically a clarity optimisation rather than a speed optimisation, and in this case I don't think it's clearer. -thomas -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. -.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. _._ ---------------------------------------------------------------------- If you have received this e-mail in error or wish to read our e-mail disclaimer statement and monitoring policy, please refer to http://www.drkw.com/disc/email/ or contact the sender. ---------------------------------------------------------------------- -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Ben Bolker wrote:> > Any better (more efficient, built-in) ideas for computing > > coef[1]+coef[2]*x+coef[3]*x^2+ ... > > than > > polynom <- function(coef,x) { > n <- length(coef) > > sum(coef*apply(matrix(c(rep(x,n),seq(0,n-1)),ncol=2),1,function(z)z[1]^z[2])) > } >You should be aware that accuracy rather than efficiency is often a problem with this representation of high order polynomials. Acton (1970) gives an example of an order 20 polynomial where extremely small changes in the coefficients of high order terms result in large changes in the roots. If your intention is to evaluate the polynomial value then you will want to consider a Horner scheme. If your intention is to find the polynomial roots then you should try to avoid the polynomial coefficient representation if at all possible. This problem comes up in ARMA and VAR time series models, where stability of the models is determined by the roots of the determinant of the AR polynomial matrix. For multivariate models this can be a fairly high order polynomial. There is a much better way to calculate the roots in this case. Paul Gilbert ----- Acton, F. S. (1970) Numerical Methods that Work, Harper & Row, NY. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._