<p>The method I wrote as of my latest edit is now faster than even <code>scipy.statstools.acf</code> with <code>fft=True</code> until the sample size gets very large.</p>
<p><strong>Error analysis</strong> If you want to adjust for biases &amp; get highly accurate error estimates: Look at my code <a href="https://github.com/flipdazed/Hybrid-Monte-Carlo/blob/master/correlations/errors.py" rel="nofollow noreferrer">here</a> which implements <a href="http://arxiv.org/pdf/hep-lat/0306017.pdf" rel="nofollow noreferrer">this paper</a> by Ulli Wolff <em>(<a href="https://www.physik.hu-berlin.de/de/com/ALPHAsoft" rel="nofollow noreferrer">or original by UW in <code>Matlab</code></a>)</em></p>
<h2>Functions Tested</h2>
<ul>
<li><code>a = correlatedData(n=10000)</code> is from a routine found <a href="https://github.com/flipdazed/Hybrid-Monte-Carlo/blob/master/correlations/acorr.py" rel="nofollow noreferrer">here</a></li>
<li><code>gamma()</code> is from same place as <code>correlated_data()</code></li>
<li><code>acorr()</code> is my function below</li>
<li><code>estimated_autocorrelation</code> is found in another answer</li>
<li><code>acf()</code> is from <code>from statsmodels.tsa.stattools import acf</code></li>
</ul>
<h2>Timings</h2>
<pre><code>%timeit a0, junk, junk = gamma(a, f=0)                            # puwr.py
%timeit a1 = [acorr(a, m, i) for i in range(l)]                   # my own
%timeit a2 = acf(a)                                               # statstools
%timeit a3 = estimated_autocorrelation(a)                         # numpy
%timeit a4 = acf(a, fft=True)                                     # stats FFT

## -- End pasted text --
100 loops, best of 3: 7.18 ms per loop
100 loops, best of 3: 2.15 ms per loop
10 loops, best of 3: 88.3 ms per loop
10 loops, best of 3: 87.6 ms per loop
100 loops, best of 3: 3.33 ms per loop
</code></pre>
<p>Edit... I checked again keeping <code>l=40</code> and changing <code>n=10000</code> to <code>n=200000</code> samples the FFT methods start to get a bit of traction and <code>statsmodels</code> fft implementation just edges it... (order is the same)</p>
<pre><code>## -- End pasted text --
10 loops, best of 3: 86.2 ms per loop
10 loops, best of 3: 69.5 ms per loop
1 loops, best of 3: 16.2 s per loop
1 loops, best of 3: 16.3 s per loop
10 loops, best of 3: 52.3 ms per loop
</code></pre>
<p>Edit 2: I changed my routine and re-tested vs. the FFT for <code>n=10000</code> and <code>n=20000</code></p>
<pre><code>a = correlatedData(n=200000); b=correlatedData(n=10000)
m = a.mean(); rng = np.arange(40); mb = b.mean()
%timeit a1 = map(lambda t:acorr(a, m, t), rng)
%timeit a1 = map(lambda t:acorr.acorr(b, mb, t), rng)
%timeit a4 = acf(a, fft=True)
%timeit a4 = acf(b, fft=True)

10 loops, best of 3: 73.3 ms per loop   # acorr below
100 loops, best of 3: 2.37 ms per loop  # acorr below
10 loops, best of 3: 79.2 ms per loop   # statstools with FFT
100 loops, best of 3: 2.69 ms per loop # statstools with FFT
</code></pre>
<h2>Implementation</h2>
<pre><code>def acorr(op_samples, mean, separation, norm = 1):
    &quot;&quot;&quot;autocorrelation of a measured operator with optional normalisation
    the autocorrelation is measured over the 0th axis

    Required Inputs
        op_samples  :: np.ndarray :: the operator samples
        mean        :: float :: the mean of the operator
        separation  :: int :: the separation between HMC steps
        norm        :: float :: the autocorrelation with separation=0
    &quot;&quot;&quot;
    return ((op_samples[:op_samples.size-separation] - mean)*(op_samples[separation:]- mean)).ravel().mean() / norm
</code></pre>
<p><strong><code>4x</code> speedup</strong> can be achieved below. You must be careful to only pass <code>op_samples=a.copy()</code> as it will modify the array <code>a</code> by <code>a-=mean</code> otherwise:</p>
<pre><code>op_samples -= mean
return (op_samples[:op_samples.size-separation]*op_samples[separation:]).ravel().mean() / norm
</code></pre>
<h2>Sanity Check</h2>
<p><a href="https://i.stack.imgur.com/3lmNf.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/3lmNf.png" alt="enter image description here" /></a></p>
<h2>Example Error Analysis</h2>
<p>This is a bit out of scope but I can't be bothered to redo the figure without the integrated autocorrelation time or integration window calculation. The autocorrelations with errors are clear in the bottom plot
<a href="https://i.stack.imgur.com/3xFQm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/3xFQm.png" alt="enter image description here" /></a></p>
Estimate Autocorrelation using Python
<p>I would like to perform Autocorrelation on the signal shown below. The time between two consecutive points is 2.5ms (or a repetition rate of 400Hz).</p>

<p><img src="https://i.stack.imgur.com/UknOq.png" alt="enter image description here"></p>

<p>This is the equation for estimating autoacrrelation that I would like to use (Taken from <a href="http://en.wikipedia.org/wiki/Autocorrelation" rel="noreferrer">http://en.wikipedia.org/wiki/Autocorrelation</a>, section Estimation):</p>

<p><img src="https://i.stack.imgur.com/CujxE.png" alt="enter image description here"></p>

<p>What is the simplest method of finding the estimated autocorrelation of my data in python? Is there something similar to <code>numpy.correlate</code> that I can use? </p>

<p>Or should I just calculate the mean and variance?</p>

<hr>

<p>Edit:</p>

<p>With help from <a href="https://stackoverflow.com/users/190597/unutbu">unutbu</a>, I have written:</p>

<pre><code>from numpy import *
import numpy as N
import pylab as P

fn = 'data.txt'
x = loadtxt(fn,unpack=True,usecols=[1])
time = loadtxt(fn,unpack=True,usecols=[0]) 

def estimated_autocorrelation(x):
    n = len(x)
    variance = x.var()
    x = x-x.mean()
    r = N.correlate(x, x, mode = 'full')[-n:]
    #assert N.allclose(r, N.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
    result = r/(variance*(N.arange(n, 0, -1)))
    return result

P.plot(time,estimated_autocorrelation(x))
P.xlabel('time (s)')
P.ylabel('autocorrelation')
P.show()
</code></pre>