Applying a half-gaussian filter to binned time series data in python

user_4770853 0
Viewed 359

I am binning some time series data, I need to apply a half-normal filter to the binned data. How can I do this in python? I've provided a toy example bellow. I need Xbinned to be smoothed with a half-gaussian filter with std of 0.25 (or what ever). I'm pretty sure the half gaussian should be facing the forward time direction.

import numpy as np

X = np.random.randint(2, size=100) #example random process

bin_size =  5

Xbinned = []

for i in range(0, len(X)+1, bin_size):
    Xbinned.append(sum(X[i:i+(bin_size-1)])/bin_size)
time-series filtering python
1 Answers

How to implement half-gaussian filtering

Scipy has a function called scipy.ndimage.gaussian_filter(). It nearly implements what we want here. Unfortunately, there's no option to use a half-gaussian instead of a gaussian. However, scipy is open-source, so we can just take the source code and modify it to be a half-gaussian.

I used this source code, and removed all of the parts that are not needed for this particular case. At the end, I had this:

import scipy.ndimage

def halfgaussian_kernel1d(sigma, radius):
    """
    Computes a 1-D Half-Gaussian convolution kernel.
    """
    sigma2 = sigma * sigma
    x = np.arange(0, radius+1)
    phi_x = np.exp(-0.5 / sigma2 * x ** 2)
    phi_x = phi_x / phi_x.sum()

    return phi_x

def halfgaussian_filter1d(input, sigma, axis=-1, output=None,
                      mode="constant", cval=0.0, truncate=4.0):
    """
    Convolves a 1-D Half-Gaussian convolution kernel.
    """
    sd = float(sigma)
    # make the radius of the filter equal to truncate standard deviations
    lw = int(truncate * sd + 0.5)
    weights = halfgaussian_kernel1d(sigma, lw)
    origin = -lw // 2
    return scipy.ndimage.convolve1d(input, weights, axis, output, mode, cval, origin)

A short summary of how this works:

  1. First, it generates a convolution kernel. It uses the formula e^(-1/2 * (x/sigma)^2) to generate the gaussian distribution. It keeps going until you're 4 standard deviations away from the center.
  2. Next, it convolves that kernel against your signal. It adjusts the kernel to start at the current timestep instead of being centered on the current timestep.

Trying this on your signal, I get a result like this:

array([0.59979879, 0.6       , 0.40006707, 0.59993293, 0.79993293,
       0.40013414, 0.20006707, 0.59986586, 0.40006707, 0.4       ,
       0.99979879, 0.00033535, 0.59979879, 0.40006707, 0.00013414,
       0.59979879, 0.20013414, 0.00006707, 0.19993293, 0.59986586])

Choice of standard deviation

If you pick a standard deviation of 0.25, that is going to have almost no effect on your signal. Here are the convolution weights it uses: [0.99966465 0.00033535]. In other words, this has less than a 0.1% effect on the signal.

I'd recommend using a larger sigma value.

Off by one error

Also, I want to point out the off-by-one error here:

for i in range(0, len(X)+1, bin_size):
    Xbinned.append(sum(X[i:i+(bin_size-1)])/bin_size)

Numpy ranges are not inclusive, so a range of i to i+(bin_size-1) actually captures 4 elements, not 5.

To fix this, you can change it to this:

for i in range(0, len(X), bin_size):
    Xbinned.append(X[i:i+bin_size].mean())

(Also, I fixed an off-by-one error in the loop specification and used a numpy shortcut for finding the mean.)

user_5301600
Related