Rank inverse weighting similar to pandas.ewm()

Question

I'm looking for a method to to get a weighted average for each row in a way where the calculated weights are halved depending on all the previous values (and that continuing for each row with growing previous values).

So an output like Out[1) should look like Out[2]:

The weights should go from 1/1 to 1/2 to 1/3 to 1/4 to 1/5... etc. depending on previous numbers of observations.

>> Out[1]:
   B
0  1
1  2
2  5
3  3

Out[2]:
   B
0  1
1  1.66
2  3.455
3  3.080


# Where row 1 is calculated as: (2*1/1 + 1*1/2)/(1/1 + 1/2) = 1.66,
# and row 2 as: (5*1/1 + 2*1/2 + 1*1/3)/(1/1 + 1/2 + 1/3) = 3.455,
# row 3 as: (3*1/1 + 5*1/2 + 2*1/3 + 1*1/4)/(1/1 + 1/2 + 1/3 + 1/4) = 3.080,
# ... and that way keeping forward depending on how much observations I have.

Would it be possible to modify the pandas.ewm() function to get a solution since it has a similar approach? I just can't get my head around on how to make this weighting for each row depending on all previous rows.

Maybe someone is able to convert that into a code working with Python?

user_4013571 · Answer

You cannot do this with pandas EWM function.

This is because you cannot represent the geometric series

S_1 = 1 + 1/2 + 1/3 + ...

in the form

S_2 = 1 + (1 - alpha) + (1 - alpha)**2 + ...

This should get you what you're after

import numpy as np


def rolling_mean_harmonic(x: np.ndarray) -> np.ndarray:
    """
    Calculates the rolling mean using harmonic series weightings
    along the first dimension of ``x``

    Example:
        >>> rolling_mean_harmonic(np.array([1, 2, 5, 5]))
        array([1.        , 1.66666667, 3.45454545, 3.08      ])
    """
    out = np.empty(x.shape, dtype=np.float64)
    h = (1/np.arange(1, len(x) + 1))
    h_cumsum = h.cumsum()
    for i in range(x.shape[0]):
        out[i] = (x[:i + 1] * h[:i + 1][::-1]).sum() / h_cumsum[i]
    return out

If speed is an issue the code above can be easily improved and parallelised with Numba njit + prange

Edit...

I tried with numba and found that parallel was actually no speed up:

import numpy as np
import numba as nb


nb.njit('float64[:](float64[:])', nogil=True, fastmath=True)
def rolling_mean_harmonic_nb_noparr(x: np.ndarray) -> np.ndarray:
    """
    Calculates the rolling mean using harmonic series weightings

    Example:
        >>> rolling_mean_harmonic(np.array([1, 2, 5, 5], dtype=float))
        array([1.        , 1.66666667, 3.45454545, 3.08      ])
    """
    out = np.empty_like(x)
    h = (1/np.arange(1, len(x) + 1))
    h_cumsum = h.cumsum()
    for i in range(x.shape[0]):
        out[i] = (x[:i + 1] * h[:i + 1][::-1]).sum() / h_cumsum[i]
    return out


nb.njit('float64[:](float64[:])', parallel=True, nogil=True, fastmath=True)
def rolling_mean_harmonic_nb(x: np.ndarray) -> np.ndarray:
    """
    Calculates the rolling mean using harmonic series weightings

    Example:
        >>> rolling_mean_harmonic(np.array([1, 2, 5, 5], dtype=float))
        array([1.        , 1.66666667, 3.45454545, 3.08      ])
    """
    out = np.empty_like(x)
    h = (1/np.arange(1, len(x) + 1))
    h_cumsum = h.cumsum()
    for i in nb.prange(x.shape[0]):
        out[i] = (x[:i + 1] * h[:i + 1][::-1]).sum() / h_cumsum[i]
    return out

Rank inverse weighting similar to pandas.ewm()

1 Answers