在不使用np.correlate的情况下计算NumPy中的自协方差函数向量

问题描述 投票:0回答:1

我正在尝试创建一个程序,使用Hannan-Rissanen算法计算ARMA(p,q)自回归移动平均随机过程的样本参数。

我遇到困难的主要步骤是计算时间序列的自协方差函数。

该程序应该采用n×1维列向量Y并计算由下式给出的k×1维列向量γ^ hat:

acvf equation image

其中Ybar是Y元素的平均值。

如何有效地计算上述总和? (显然for循环可以工作,但我正试图在矢量化numpy操作上做得更好)因为我正在使用它作为一种学习经验,我宁愿不使用除了像np.sumnp.mean这样的非常基本的函数之外的任何numpy函数。

以下类似的问题已被提出,但不能完全回答我的问题:

Computing autocorrelation of vectors with numpy(使用np.correlate

(其他一些人遇到了使用更高级的numpy函数的相同问题,或者没有像我希望的那样吐出向量。)

arrays numpy
1个回答
0
投票

这是替换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]

© www.soinside.com 2019 - 2024. All rights reserved.