1

I was trying to take a oscillation avarage of a highly oscillating data. The oscillations are not uniform, it has less oscillations in the initial regions.

x = np.linspace(0, 1000, 1000001)
y = some oscillating data say, sin(x^2)

(The original data file is huge, so I can't upload it)

I want to take a weighted moving avarage of the function and plot it. Initially the period of the function is larger, so I want to take avarage over a large time interval. While I can do with smaller time interval latter.

I have found a possible elegant solution in following post:

Weighted moving average in python

However, I want to have different width in different regions of x. Say when x is between (0,100) I want the width=0.6, while when x is between (101, 300) width=0.2 and so on.

This is what I have tried to implement( with my limited knowledge in programing!)

def weighted_moving_average(x,y,step_size=0.05):#change the width to control average
    bin_centers  = np.arange(np.min(x),np.max(x)-0.5*step_size,step_size)+0.5*step_size
    bin_avg = np.zeros(len(bin_centers))

    #We're going to weight with a Gaussian function
    def gaussian(x,amp=1,mean=0,sigma=1):
         return amp*np.exp(-(x-mean)**2/(2*sigma**2))


    if x.any()  < 100:
        for index in range(0,len(bin_centers)):
            bin_center = bin_centers[index]
            weights = gaussian(x,mean=bin_center,sigma=0.6)
            bin_avg[index] = np.average(y,weights=weights)

    else:
        for index in range(0,len(bin_centers)):
            bin_center = bin_centers[index]
            weights = gaussian(x,mean=bin_center,sigma=0.1)
            bin_avg[index] = np.average(y,weights=weights)

    return (bin_centers,bin_avg)

It is needless to say that this is not working! I am getting the plot with the first value of sigma. Please help...

Snow bunting
  • 1,120
  • 8
  • 28
Archimedes
  • 365
  • 3
  • 15
  • Your question is a bit too general and math oriented for stack overflow. Try to edit your post with an explicit expected outcome (like "I want to generate a list composed from bla bla bla"). Also, the scope of your function is not clear. Is "return (bin_centers,bin_avg)" supposed to be in the function? – zar3bski Jun 09 '18 at 07:45
  • I have edited the question. You may also refer to the linked question for initial details. – Archimedes Jun 09 '18 at 08:02

1 Answers1

2

The following snippet should do more or less what you tried to do. You have mainly a logical problem in your code, x.any() < 100 will always be True, so you'll never execute the second part.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 10, 1000)
y = np.sin(x**2)

def gaussian(x,amp=1,mean=0,sigma=1):
    return amp*np.exp(-(x-mean)**2/(2*sigma**2))

def weighted_average(x,y,step_size=0.3):
    weights = np.zeros_like(x)
    bin_centers  = np.arange(np.min(x),np.max(x)-.5*step_size,step_size)+.5*step_size
    bin_avg = np.zeros_like(bin_centers)

    for i, center in enumerate(bin_centers):
        # Select the indices that should count to that bin
        idx = ((x >= center-.5*step_size) & (x <= center+.5*step_size))

        weights = gaussian(x[idx], mean=center, sigma=step_size)
        bin_avg[i] = np.average(y[idx], weights=weights)
    return (bin_centers,bin_avg)

idx = x <= 4
plt.plot(*weighted_average(x[idx],y[idx], step_size=0.6))

idx = x >= 3
plt.plot(*weighted_average(x[idx],y[idx], step_size=0.1))

plt.plot(x,y)
plt.legend(['0.6', '0.1', 'y'])
plt.show()

enter image description here

However, depending on the usage, you could also implement moving average directly:

x = np.linspace(0, 60, 1000)
y = np.sin(x**2)

z = np.zeros_like(x)
z[0] = x[0]

for i, t in enumerate(x[1:]):
    a=.2
    z[i+1] = a*y[i+1] + (1-a)*z[i]

plt.plot(x,y)
plt.plot(x,z)
plt.legend(['data', 'moving average'])
plt.show()

enter image description here

Of course you could then change a adaptively, e.g. depending of the local variance. Also note that this has apriori a small bias depending on a and the step size in x.

Snow bunting
  • 1,120
  • 8
  • 28