Оценка автокорреляции с использованием Python

Ответ 1

Я не думаю, что для этого конкретного вычисления есть функция NumPy. Вот как я его написал:

def estimated_autocorrelation(x):
    """
    http://stackoverflow.com/q/14297012/190597
    http://en.wikipedia.org/wiki/Autocorrelation#Estimation
    """
    n = len(x)
    variance = x.var()
    x = x-x.mean()
    r = np.correlate(x, x, mode = 'full')[-n:]
    assert np.allclose(r, np.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
    result = r/(variance*(np.arange(n, 0, -1)))
    return result

Утверждающий оператор существует как для проверки расчета, так и для документирования его намерения.

Если вы уверены, что эта функция ведет себя как ожидалось, вы можете прокомментировать инструкцию assert или запустить script с помощью python -O. (Флаг -O указывает Python игнорировать утверждения assert.)

Ответ 2

Я взял часть кода из функции pandas autocorrelation_plot(). Я проверил ответы с помощью R, и значения точно совпадают.

import numpy
def acf(series):
    n = len(series)
    data = numpy.asarray(series)
    mean = numpy.mean(data)
    c0 = numpy.sum((data - mean) ** 2) / float(n)

    def r(h):
        acf_lag = ((data[:n - h] - mean) * (data[h:] - mean)).sum() / float(n) / c0
        return round(acf_lag, 3)
    x = numpy.arange(n) # Avoiding lag 0 calculation
    acf_coeffs = map(r, x)
    return acf_coeffs

Ответ 3

Пакет statsmodels добавляет функцию автокорреляции, которая внутренне использует np.correlate (согласно документации statsmodels).

См: http://statsmodels.sourceforge.net/stable/generated/statsmodels.tsa.stattools.acf.html#statsmodels.tsa.stattools.acf

Ответ 4

Метод, который я написал с моего последнего редактирования, теперь быстрее, чем даже scipy.statstools.acf с fft=True, пока размер выборки не станет очень большим.

Анализ ошибок. Если вы хотите настроить отклонения и получить очень точные оценки ошибок: посмотрите на мой код здесь который реализует эту статью Улли Вольфа (или оригинал UW в Matlab)

Проверенные функции

• a = correlatedData(n=10000) происходит из найденной процедуры здесь
• gamma() - это то же место, что и correlated_data()
• acorr() - моя функция ниже
• estimated_autocorrelation находится в другом ответе
• acf() - от from statsmodels.tsa.stattools import acf

Задержки

%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

Изменить... Я снова проверил, сохраняя l=40 и меняя n=10000 на n=200000 образцы, методы FFT начинают получать немного тяги, а реализация statsmodels fft просто режет его... (заказ - это то же самое)

## -- 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

Редактировать 2: я изменил свою процедуру и повторно протестировал против FFT для n=10000 и n=20000

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

Реализация

def acorr(op_samples, mean, separation, norm = 1):
    """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
    """
    return ((op_samples[:op_samples.size-separation] - mean)*(op_samples[separation:]- mean)).ravel().mean() / norm

4x ускорение может быть достигнуто ниже. Вы должны быть осторожны, чтобы пройти только op_samples=a.copy(), так как он изменит массив a на a-=mean иначе:

op_samples -= mean
return (op_samples[:op_samples.size-separation]*op_samples[separation:]).ravel().mean() / norm

Проверка работоспособности

Пример анализа ошибок

Это немного выходит за пределы области видимости, но я не могу потрудиться, чтобы переделать фигуру без встроенного расчета автокорреляции или расчета окна интеграции. Автокорреляции с ошибками ясны в нижнем графике

Ответ 5

Я нашел, что полученные ожидаемые результаты произошли с небольшим изменением:

def estimated_autocorrelation(x):
    n = len(x)
    variance = x.var()
    x = x-x.mean()
    r = N.correlate(x, x, mode = 'full')
    result = r/(variance*n)
    return result

Тестирование результатов автокорреляции Excel.