Can the cumsum function in NumPy decay while adding?

Question

I have an array of values a = (2,3,0,0,4,3)

y=0
for x in a:
  y = (y+x)*.95

Is there any way to use cumsum in numpy and apply the .95 decay to each row before adding the next value?

Cass

Answer 1

You're asking for a simple IIR Filter. Scipy's lfilter() is made for that:

import numpy as np
from scipy.signal import lfilter

data = np.array([2, 3, 0, 0, 4, 3], dtype=float)  # lfilter wants floats

# Conventional approach:
result_conv = []
last_value = 0
for elmt in data:
    last_value = (last_value + elmt)*.95
    result_conv.append(last_value)

# IIR Filter:
result_IIR = lfilter([.95], [1, -.95], data)

if np.allclose(result_IIR, result_conv, 1e-12):
    print("Values are equal.")

Disaffection

I'd consider using xx = np.array([2, 3, 0, 0, 4, 3], dtype=np.float64) instead here – Gove

Valid point. I reran my code and it seems that the complaint from lfilter() came from a different problem. – Disaffection

Answer 2

If you're only dealing with a 1D array, then short of scipy conveniences or writing a custom reduce ufunc for numpy, then in Python 3.3+, you can use itertools.accumulate, eg:

from itertools import accumulate

a = (2,3,0,0,4,3)
y = list(accumulate(a, lambda x,y: (x+y)*0.95))
# [2, 4.75, 4.5125, 4.286875, 7.87253125, 10.3289046875]

Gove

Answer 3

Numba provides an easy way to vectorize a function, creating a universal function (thus providing ufunc.accumulate):

import numpy
from numba import vectorize, float64

@vectorize([float64(float64, float64)])
def f(x, y):
    return 0.95 * (x + y)

>>> a = numpy.array([2, 3, 0, 0, 4, 3])
>>> f.accumulate(a)
array([  2.        ,   4.75      ,   4.5125    ,   4.286875  ,
         7.87253125,  10.32890469])

Deodorize

Answer 4

I don't think that this can be done easily in NumPy alone, without using a loop.

One array-based idea would be to calculate the matrix M_ij = .95**i * a[N-j] (where N is the number of elements in a). The numbers that you are looking for are found by summing entries diagonally (with i-j constant). You could use thus use multiple numpy.diagonal(…).sum().

The good old algorithm that you outline is clearer and probably quite fast already (otherwise you can use Cython).

Doing what you want through NumPy without a single loop sounds like wizardry to me. Hats off to anybody who can pull this off.

Martinet