Вычисление вектора автоковариантной функции в NumPy без использования np.correlate

Question

Я пытаюсь создать программу, использующую алгоритм Ханнана-Риссанена для вычисления параметров выборки для авторегрессионного процесса скользящей средней ARMA (p, q).

Основной шаг, с которым я сталкиваюсь, - это вычисление автоковариационной функции временного ряда.

Программа должна взять n × 1-мерный вектор столбца Y и вычислить ak × 1.-мерный вектор-столбец γ ^ hat, определяемый как:

изображение уравнения acvf

где Ybar - среднее значение элементов Y.

Как можноЯ рассчитываю вышеуказанную сумму эффективно?(очевидно, что циклы будут работать, но я пытаюсь лучше справляться с векторизованными операциями с NumPy). Поскольку я использую это в качестве обучающего опыта, я бы предпочел не использовать какие-либо функции NUMPY, кроме очень простых, таких как np.sum илиnp.mean.

Был задан следующий предыдущий аналогичный вопрос, но он не вполне отвечает на мой вопрос:

Вычисление автокорреляции векторов с numpy (используется * 1021)*)

(некоторые другие страдают от той же проблемы использования более продвинутых numpy функций или не выделяют векторы, как я хотел бы сделать здесь.)

Answer 1

Вот один из способов замены np.correlate (я полагаю, что это основная трудность; я также предполагаю, что у вас нет желания сдавать код в FFT):

def autocorr_direct(a, debug=False):
    n, _ = a.shape
    out = np.zeros((n+1, 2*n-1), a.dtype)
    out.reshape(-1)[:2*n*n].reshape(n, 2*n)[::-1, :n] = a*a.T
    if debug:
        print(out.reshape(-1)[:2*n*n].reshape(n, 2*n))
        print(out)
    return out.sum(0)

Например:

>>> a = np.array([[1, 1, 2, -1]]).T
>>> autocorr_direct(a, True)
[[-1 -1 -2  1  0  0  0  0]
 [ 2  2  4 -2  0  0  0  0]
 [ 1  1  2 -1  0  0  0  0]
 [ 1  1  2 -1  0  0  0  0]]
[[-1 -1 -2  1  0  0  0]
 [ 0  2  2  4 -2  0  0]
 [ 0  0  1  1  2 -1  0]
 [ 0  0  0  1  1  2 -1]
 [ 0  0  0  0  0  0  0]]
array([-1,  1,  1,  7,  1,  1, -1])
>>> np.correlate(a[:, 0], a[:, 0], 'full')
array([-1,  1,  1,  7,  1,  1, -1])

Обратите внимание на трюк с изменением формы, который срезает квадратный массив a[::-1]*a.T.

Примечание 2; чтобы получить вектор столбца из 1D вектора X используйте X[:, None].