[R] polynomial

[Paul Gilbert]

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