[Numpy-discussion] poly class question
Sebastian Walter
sebastian.walter@gmail....
Mon Oct 5 04:37:39 CDT 2009
On Fri, Oct 2, 2009 at 10:40 PM, <josef.pktd@gmail.com> wrote:
> On Fri, Oct 2, 2009 at 3:38 PM, Charles R Harris
> <charlesr.harris@gmail.com> wrote:
>>
>>
>> On Fri, Oct 2, 2009 at 12:30 PM, <josef.pktd@gmail.com> wrote:
>>>
>>> On Fri, Oct 2, 2009 at 2:09 PM, Charles R Harris
>>> <charlesr.harris@gmail.com> wrote:
>>> >
>>> >
>>> > On Fri, Oct 2, 2009 at 11:35 AM, Charles R Harris
>>> > <charlesr.harris@gmail.com> wrote:
>>> >>
>>> >>
>>> >> On Fri, Oct 2, 2009 at 11:33 AM, Charles R Harris
>>> >> <charlesr.harris@gmail.com> wrote:
>>> >>>
>>> >>>
>>> >>> On Fri, Oct 2, 2009 at 11:30 AM, Charles R Harris
>>> >>> <charlesr.harris@gmail.com> wrote:
>>> >>>>
>>> >>>>
>>> >>>> On Fri, Oct 2, 2009 at 11:08 AM, <josef.pktd@gmail.com> wrote:
>>> >>>>>
>>> >>>>> Is there a way in numpy (or scipy) to get an infinite expansion for
>>> >>>>> the inverse of a polynomial (for a finite number of terms)
>>> >>>>>
>>> >>>>> np.poly1d([ -0.8, 1])**(-1)
>>> >>>>>
>>> >>>>> application for example the MA representation of an AR(1)
>>> >>>>>
>>> >>>>
>>> >>>> Hmm, I've been working on a chebyshev class and division of a scalar
>>> >>>> by
>>> >>>> a chebyshev series is
>>> >>>> expressly forbidden, but it could be included if a good interface is
>>> >>>> proposed. Same would go for polynomials.
>>> >>>
>>> >>> In fact is isn't hard to get, for poly1d you should be able to
>>> >>> multiply
>>> >>> the series by a power of x to shift it left, then divide.
>>> >>>
>>> >>
>>> >> That is, divide a power of x by the polynomial.
>>> >>
>>> >
>>> > You will also need to reverse the denominator coefficients...Chuck
>>>
>>> That's the hint I needed. However the polynomial coefficients are then
>>> reversed and not consistent with other polynomial operations, aren't
>>> they?
>>>
>>> >>> from scipy.signal import lfilter
>>>
>>> >>> (np.poly1d([1, 0])**10)/np.poly1d([1, -0.8])
>>> (poly1d([ 1. , 0.8 , 0.64 , 0.512 , 0.4096 ,
>>> 0.32768 , 0.262144 , 0.2097152 , 0.16777216,
>>> 0.13421773]), poly1d([ 0.10737418]))
>>>
>>> >>> lfilter([1], [1,-0.8], [1] + [0]*9)
>>> array([ 1. , 0.8 , 0.64 , 0.512 , 0.4096 ,
>>> 0.32768 , 0.262144 , 0.2097152 , 0.16777216, 0.13421773])
>>>
>>> >>> (np.poly1d([1, 0])**10)/np.poly1d([1, -0.8, 0.2])
>>> (poly1d([ 1. , 0.8 , 0.44 , 0.192 , 0.0656 ,
>>> 0.01408 , -0.001856 , -0.0043008 , -0.00306944]),
>>> poly1d([-0.00159539, 0.00061389]))
>>> >>> lfilter([1], [1,-0.8, 0.2], [1] + [0]*9)
>>> array([ 1. , 0.8 , 0.44 , 0.192 , 0.0656 ,
>>> 0.01408 , -0.001856 , -0.0043008 , -0.00306944, -0.00159539])
>>>
>>>
>>> What I meant initally doesn't necessarily mean division of a scalar.
>>>
>>> >>> np.poly1d([1])/np.poly1d([-0.8, 1])
>>> (poly1d([ 0.]), poly1d([ 1.]))
>>>
>>> I didn't find any polynomial division that does the expansion of the
>>> remainder. The same problem, I think is inherited, by the
>>> scipy.signal.lti, and it took me a while to find the usefulness of
>>> lfilter in this case.
>>>
>>> If it were possible to extend the methods for the polynomial class to
>>> do a longer expansions, it would make them more useful for arma and
>>> lti.
>>>
>>> (in some areas, I'm still trying to figure out whether some
>>> functionality is just hidden to me, or actually a limitation of the
>>> implementation or a missing feature.)
>>>
>>
>> Could you describe the sort of problems you want to solve? There are lots of
>> curious things out there we could maybe work with. Covariances, for
>> instance, are closely related to Chebyshev series.
>
> I am working on a discrete time arma process of the form
>
> a(L) x_t = b(L) u_t, where L is the lag operator L^k x_t = x_(t-k)
>
> what I just programmed using lfilter is
> x_t = b(L)/a(L) u_t where b(L)/a(L) is the impulse response function
> or moving average representation
>
> a(L)/b(L) is the autoregressive representation
>
> the extension
> a(L)(1-L)^d x_t = b(L) u_t, where d = 0,1,2,... (standard) or also
> continuous d <1 (fractional integration)
>
> a(L)/b(L), b(L)/a(L) (1-L)^(-d) or (1-L)^d (0<d<1) are infinite
> dimensional lag polynomials in the general case.
> Initially I was looking for an easy way to do these calculation as polynomials.
> (The fractional case (1-L)^d (0<d<1) might be a pretty special case,
> and I just looked it up today, but is a popular model class in
> econometrics, fractionally integrated arma processes)
>
> multiplication works well np.poly1d([-1, 1])*np.poly1d([-0.8, 1])
> (with reversed poly coefficient scipy.signal I think)
>
> the functions in scipy.signal for lti are only for continuous time
> processes and use poly1d under the hood, which means for example
>
>>>> from scipy import signal
>>>> signal.impulse(([1, -0.8],[1]), N=10)
> raise ValueError, "Improper transfer function."
> ValueError: Improper transfer function.
>
> while this works
>>>> signal.impulse(([1],[1, -0.8]), N=10)
>
> (It's been a while since I looked inside scipy.signal.lti)
>
> A separate issue would be the multivariate version VARMA, or MIMO in
> system modeling. a(L), b(L) are matrix polynomials and x_t, u_t are
> 1d arrays evolving in time.
> But that is a different discussion.
I'm working on something that requires truncated
univariate/multivariate operations
on scalar, vector and matrix polynomials. I'm pretty sure I'm doing
something different
than what is used in MIMO. Still, do you have a good reference for
implementation/complexity/stability of operations on multivariate
(matrix) polynomials?
Want to make sure I'm not reinventing the wheel ;)
>
> I'm not very familiar with Chebychev polynomials, the last time I
> wanted to use them I didn't see anything about their use as a base for
> functions in several variables and gave up. I've seen papers that use
> them as base for functions in one variable, but I'm not doing anything
> like this right now.
>
> Thanks,
>
> Josef
>
>>
>> Chuck
>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
More information about the NumPy-Discussion
mailing list