4
\$\begingroup\$

I'm implementing the Probabilistic Exponentially Weighted Mean for real time prediction of sensor data in pandas but have issues with optimising the pandas notebook for quick iterations.

Is there a more optimal way to completely remove the for loop as it currently runs longer than expected. How can I take advantage of apply() ; or vectorized operations etc. here?

import pandas as pd
import numpy as np
from datetime import datetime, timedelta
import datetime
import matplotlib.pyplot as plt
#plt.style.use('fivethirtyeight')
#%config InlineBackend.figure_format = 'retina'
#%matplotlib inline
from itertools import islice
from math import sqrt
from scipy.stats import norm


ts = pd.date_range(start ='1-1-2019',  
         end ='1-10-2019', freq ='5T') 


np.random.seed(seed=1111)
data = np.random.normal(2.012547, 1.557331,size=len(ts))
df = pd.DataFrame({'timestamp': ts, 'speed': data})
df.speed = df.speed.abs()
df = df.set_index('timestamp')
time_col = 'timestamp'
value_col = 'speed'

#pewna parameters
T = 30      # initialization period (in cycles)
beta = 0.5  # lower values make the algorithm behave more like regular EWMA
a = 0.99    # the maximum value of the EWMA a parameter, used for outliers
z = 3


#the PEWNA Model

# create a DataFrame for the run time variables we'll need to calculate
pewm = pd.DataFrame(index=df.index, columns=['Mean', 'Var', 'Std'], dtype=float)
pewm.iloc[0] = [df.iloc[0][value_col], 0, 0]

t = 0

for _, row in islice(df.iterrows(), 1, None):
    diff = row[value_col] - pewm.iloc[t].Mean # difference from moving average
    p = norm.pdf(diff / pewm.iloc[t].Std) if pewm.iloc[t].Std != 0 else 0 # Prob of observing diff
    a_t = a * (1 - beta * p) if t > T else 1 - 1/(t+1) # weight to give to this point
    incr = (1 - a_t) * diff

    # Update Mean, Var, Std
    pewm.iloc[t+1].Mean = pewm.iloc[t].Mean + incr
    pewm.iloc[t+1].Var = a_t * (pewm.iloc[t].Var + diff * incr)
    pewm.iloc[t+1].Std = sqrt(pewm.iloc[t+1].Var)
    t += 1

\$\endgroup\$
  • 1
    \$\begingroup\$ @Graipher Ive added dummy data to show what the df looks like. \$\endgroup\$ – Matimba Jun 10 '19 at 9:24
2
\$\begingroup\$

I think you can use just raw numpy arrays for speeding it up:

after #the PEWMA Model you can just use following code:

_x = df[value_col]
_mean, _std, _var = np.zeros(_x.shape), np.zeros(_x.shape), np.zeros(_x.shape)

for i in range(1, len(_x)):
    diff = _x[i] - _mean[i-1]

    p = norm.pdf(diff / _std[i-1]) if _std[i-1] != 0 else 0 # Prob of observing diff
    a_t = a * (1 - beta * p) if (i-1) > T else 1 - 1/i # weight to give to this point
    incr = (1 - a_t) * diff

    # Update Mean, Var, Std
    v = a_t * (_var[i-1] + diff * incr)
    _mean[i] = _mean[i-1] + incr
    _var[i] = v
    _std[i] = np.sqrt(v)

pewm = pd.DataFrame({'Mean': _mean, 'Var': _var, 'Std': _std}, index=df.index)

If you need to apply it to really big data you could make it even faster by using numba package:

from numba import njit

# numba has some issues with stats's norm.pdf so redefine it here as function
@njit
def norm_pdf(x):
    return np.exp(-x**2/2)/np.sqrt(2*np.pi)

@njit
def pwma(_x, a, beta, T):
    _mean, _std, _var = np.zeros(_x.shape), np.zeros(_x.shape), np.zeros(_x.shape)
    _mean[0] = _x[0]

    for i in range(1, len(_x)):
        diff = _x[i] - _mean[i-1]

        p = norm_pdf(diff / _std[i-1]) if _std[i-1] != 0 else 0 # Prob of observing diff
        a_t = a * (1 - beta * p) if (i-1) > T else 1 - 1/(i) # weight to give to this point
        incr = (1 - a_t) * diff

        # Update Mean, Var, Std
        v = a_t * (_var[i-1] + diff * incr)
        _mean[i] = _mean[i-1] + incr
        _var[i] = v
        _std[i] = np.sqrt(v)
    return _mean, _var, _std

# Using :
_mean, _var, _std = pwma(df[value_col].values, 0.99, 0.5, 30)
pewm = pd.DataFrame({'Mean': _mean, 'Var': _var, 'Std': _std}, index=df.index)
| improve this answer | |
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.