CODEWITHSUNDEEP
pythonarraysnumpycumsum
D_C_A
D_C_A

Reputation: 111

Can the cumsum function in NumPy decay while adding?

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?

Upvotes: 11

Views: 1635

Answers (4)

Josh Bode
Josh Bode

Reputation: 3742

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

Upvotes: 3

Dietrich
Dietrich

Reputation: 5531

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

Upvotes: 8

Eric O. Lebigot
Eric O. Lebigot

Reputation: 94565

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.

Upvotes: 1

Jon Clements
Jon Clements

Reputation: 142206

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]

Upvotes: 5

Related Questions