[Numpy-discussion] polynomial fromroots

[Charles R Harris]

On Sat, Oct 9, 2010 at 8:36 PM, <josef.pktd at gmail.com> wrote:

> On Sat, Oct 9, 2010 at 10:01 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
> >
> >
> > On Sat, Oct 9, 2010 at 7:47 PM, <josef.pktd at gmail.com> wrote:
> >>
> >> I'm trying to see whether I can do this without reading the full manual.
> >>
> >> Is it intended that fromroots normalizes the highest order term
> >> instead of the lowest?
> >>
> >>
> >> >>> import numpy.polynomial as poly
> >>
> >> >>> p = poly.Polynomial([1, -1.88494037,  0.0178126 ])
> >> >>> p
> >> Polynomial([ 1.        , -1.88494037,  0.0178126 ], [-1.,  1.])
> >> >>> pr = p.roots()
> >> >>> pr
> >> array([   0.53320748,  105.28741219])
> >> >>> poly.Polynomial.fromroots(pr)
> >> Polynomial([  56.14003571, -105.82061967,    1.        ], [-1.,  1.])
> >> >>>
> >>
> >> renormalizing
> >>
> >> >>> p2 = poly.Polynomial.fromroots(pr)
> >> >>> p2/p2.coef[0]
> >> Polynomial([ 1.        , -1.88494037,  0.0178126 ], [-1.,  1.])
> >>
> >>
> >> this is, I think what I want to do, invert roots that are
> >> inside/outside the unit circle (whatever that means
> >>
> >> >>> pr[np.abs(pr)<1] = 1./pr[np.abs(pr)<1]
> >> >>> p3 = poly.Polynomial.fromroots(pr)
> >> >>> p3/p3.coef[0]
> >> Polynomial([ 1.        , -0.54270529,  0.0050643 ], [-1.,  1.])
> >>
> >
> > Wrong function ;) You defined the polynomial by its coefficients. What
> you
> > want to do is
>
> My coefficients are from a lag-polynomial in time series analysis
> (ARMA), and they really are the (estimated) coefficients. It is
> essentially the same as the model for scipy.signal.lfilter.
> I just need to check the roots to see whether the process is
> stationary and invertible. If one of the two lag-polynomials (moving
> average) has roots on the wrong side of the unit circle, then I can
> invert them.
>
> I'm coding from memory of how this is supposed to work, so maybe I'm
> back to RTFM and RTFTB (TB=text book).
>
> (I think what I really would need is a z-transform, but I don't have
> much of an idea how to do this on a computer)
>
> Thanks, the main thing I need to do is check the convention or
> definition for the normalization. And as btw, I like that the coef are
> in increasing order
> e.g. seasonal differencing multiplied with 1 lag autoregressive
> poly.Polynomial([1.,0,0,-1])*poly.Polynomial([1,0.8])
>
>
The polynomial is obtained from the product of the terms (x - r_i), so it is
monic. In a different basis it may not appear that way:

In [23]: p = poly.Legendre.fromroots([1, -1.88494037,  0.0178126 ])

In [24]: p
Out[24]: Legendre([ 0.32261828, -1.30070346,  0.57808518,  0.4       ],
[-1.,  1.])

The z transform of a sequence is simply the polynomial with 1/z in place of
x, so you just need to invert all the roots. However, the easy way to look
at it is that you are finding  \alpha such that  \alpha^n is a solution to
the homogeneous equation. For instance, if y_n = a*y_{n - 1} + b*y_{n-2},
and you just substitute and divide out the extraneous powers, you will get
\alpha^2 = a*\alpha + b. If one of the roots has a modulus bigger than 1
then it will grow geometrically with n, meaning that roundoff errors are
likely to overwhelm the result. That is the unstable case. It is actually
the discrete version of solving homogeneous linear differential equations
with constant coefficients, so repeated roots need some care.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20101009/cf506a88/attachment.html>