kanayamalakar kanayamalakar - 7 months ago 214
Python Question

Calculating Autocorrelation function with python

All that I am trying to do is to calculate the auto correlation of an array jx for which I am using the following formula,

autocorrelation formula

where n is the time at which I wish to calculate the autocorrelation function,

Mt
is the maximum time and
tk
are time steps running from
1
to
Mt-n
.

This is the code that I have written. I am checking my program with a simple array
jx=linspace(1,10,20)
. I am also making the program save the autocorrelation values for different
n
and plot them with
n
.

from numpy import *
from pylab import*

jx=linspace(1,10,20)

Mt=len(jx)

def Hcacf(n):
Sum=0.0
coeff1=0
while coeff1 < (Mt-n) :
Sum = Sum + jx[coeff1]*jx[coeff1+n]# + jy[coeff1]*jy[coeff1+n]
coeff1=coeff1+1
avg = Sum*1.0 / (Mt-n)
return avg

autocorrelation=[]
for n in linspace(0,Mt-1,Mt):
ac=Hcacf(n+1)
autocorrelation.append(ac)

lag=linspace(0,Mt-1,Mt)
plot(lag,autocorrelation,marker='o')
show()


The output it returns is this:
enter image description here

But it also returns the following error messsage:

RuntimeWarning: invalid value encountered in double_scalars
avg = Sum*1.0 / (Mt-n)


Where am I going wrong?

Answer

It seems that you have a division by zero. It works this way:

1) In the line for n in linspace(0,Mt-1,Mt): you have n==Mt-1,

2) So, in the next line ac=Hcacf(n+1) you call the function Hcacf(Mt),

3) But inside this function Hcacf you have a line avg = Sum*1.0 / (Mt-n). It is the place where division by zero is possible.

To fix it you can exclude the endpoint of the interval in the first line. Try to replace it by this line:

for n in linspace(0, Mt-1, num=Mt, endpoint=False):