从时间序列数据计算转移矩阵的有效方法是什么？

Question

我正在尝试从时间序列数据中计算转换矩阵。我编写了一个自定义函数，如下面的代码，用于我的目的。

def compute_transition_matrix(data, n, step = 1):
    P = np.zeros((n, n))
    m = len(data)
    for i in range(m):
        initial, final = i, i + step
        if final < m:
            P[data[initial]][data[final]] += 1
    sums = np.sum(P, axis = 1)
    for i in range(n):
        for j in range(n):
            P[i][j] = P[i][j] / sums[i]
    return P

print(compute_transition_matrix([3, 0, 1, 3, 2, 6, 5, 4, 7, 5, 4], 8, 1))

上述函数中，data为输入时序数据，n为马尔可夫链中状态总数，step为转移步。

作为一个示例，我采取了，

data = [3, 0, 1, 3, 2, 6, 5, 4, 7, 5, 4]
n = 8 (this means there are 8 states in Markov chain from 0 - 7, both inclusive)
step = 1

但是，我只是想知道是否有办法使用 NumPy/pandas/scikit 中的内置函数来实现这一点？

Answer 1

我不确定是否有内置函数可以实现这一点，但我可以考虑在numpy（使用fancy indexing,broadcasting 和stride tricks）中这样做：

def compute_transition_matrix2(data, n, step = 1):
    
    t = np.array(data)
    step = step
    total_inds = t.size - (step + 1) + 1
    t_strided = np.lib.stride_tricks.as_strided(
                                    t,
                                    shape = (total_inds, 2),
                                    strides = (t.strides[0], step * t.strides[0]))
    
    inds, counts = np.unique(t_strided, axis = 0, return_counts = True)

    P = np.zeros((n, n))
    P[inds[:, 0], inds[:, 1]] = counts
    
    sums = P.sum(axis = 1)
    # Avoid divide by zero error by normalizing only non-zero rows
    P[sums != 0] = P[sums != 0] / sums[sums != 0][:, None]
    
    # P = P / P.sum(axis = 1)[:, None]
    return P

print(compute_transition_matrix2([3, 0, 1, 3, 2, 6, 5, 4, 7, 5, 4], 8, 1))

[[0.  1.  0.  0.  0.  0.  0.  0. ]
 [0.  0.  0.  1.  0.  0.  0.  0. ]
 [0.  0.  0.  0.  0.  0.  1.  0. ]
 [0.5 0.  0.5 0.  0.  0.  0.  0. ]
 [0.  0.  0.  0.  0.  0.  0.  1. ]
 [0.  0.  0.  0.  1.  0.  0.  0. ]
 [0.  0.  0.  0.  0.  1.  0.  0. ]
 [0.  0.  0.  0.  0.  1.  0.  0. ]]

您的代码结果：

def compute_transition_matrix(data, n, step = 1):
    P = np.zeros((n, n))
    m = len(data)
    for i in range(m):
        initial, final = i, i + step
        if final < m:
            P[data[initial]][data[final]] += 1
    sums = np.sum(P, axis = 1)
    for i in range(n):
        if sums[i] != 0: # Added this check
            for j in range(n):
                P[i][j] = P[i][j] / sums[i]
    return P

print(compute_transition_matrix([3, 0, 1, 3, 2, 6, 5, 4, 7, 5, 4], 8, 1))

[[0.  1.  0.  0.  0.  0.  0.  0. ]
 [0.  0.  0.  1.  0.  0.  0.  0. ]
 [0.  0.  0.  0.  0.  0.  1.  0. ]
 [0.5 0.  0.5 0.  0.  0.  0.  0. ]
 [0.  0.  0.  0.  0.  0.  0.  1. ]
 [0.  0.  0.  0.  1.  0.  0.  0. ]
 [0.  0.  0.  0.  0.  1.  0.  0. ]
 [0.  0.  0.  0.  0.  1.  0.  0. ]]

我的代码中的中间值：（供您参考）

t_strided =

array([[3, 0],
       [0, 1],
       [1, 3],
       [3, 2],
       [2, 6],
       [6, 5],
       [5, 4],
       [4, 7],
       [7, 5],
       [5, 4]])

inds, counts =

(array([[0, 1],
        [1, 3],
        [2, 6],
        [3, 0],
        [3, 2],
        [4, 7],
        [5, 4],
        [6, 5],
        [7, 5]]),
 array([1, 1, 1, 1, 1, 1, 2, 1, 1]))

时间比较：

# Generate some random large data
n = 1000
t = np.random.choice(np.arange(n), size = n)
data = list(t)

%timeit compute_transition_matrix(data, n, 1)
# 433 ms ± 21.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit compute_transition_matrix2(data, n, 1)
# 5.5 ms ± 304 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)