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?

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.")

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]

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])

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.