计算1和0序列(马尔可夫链)的智能方法

Question

假设1和0的序列(np.array)。我需要构造二阶马尔可夫链。为此，我需要计算T001的比例，…T111，其中T001是序列出现[0,0,1]的次数。显然我可以做一个for循环，但我尽量避免这样做。理想情况下，该方法可以扩展到k阶。

Answer 1

该任务可以分为两个子任务:第一个是创建长度为3的所有子序列，然后计算所有唯一子序列。这两个任务都可以通过内置的numpy方法完成:

假设输入序列在a中存储为numpy数组，则第一点可以通过

实现

windows = np.lib.stride_tricks.sliding_window_view(a, 3)

然后我们可以使用 来计算不同的窗口。

unique_windows, counts = np.unique(windows, axis=0, return_counts=True)

Answer 2

您可以使用more_itertools.windowed和collections.Counter的组合:

from more_itertools import windowed
from collections import Counter
import numpy as np
np.random.seed(0)
a = np.random.choice([0,1], size=100)
Counter(windowed(a, n=3))

输入:

array([0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1,
1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0,
0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1,
1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0])

输出:

Counter({(0, 1, 1): 14,
(1, 1, 0): 14,
(1, 0, 1): 15,
(1, 1, 1): 16,
(1, 0, 0): 10,
(0, 0, 1): 10,
(0, 1, 0): 12,
(0, 0, 0): 7})

带自定义键:

{'T'+''.join(map(str, k)): v for k,v in Counter(windowed(a, n=3)).items()}

输出:

{'T011': 14,
'T110': 14,
'T101': 15,
'T111': 16,
'T100': 10,
'T001': 10,
'T010': 12,
'T000': 7}

如何确保所有值都存在(即使计数为零):

简单的解决方案:生成一个包含所有可能键的列表，并使用它来迭代

from itertools import product
possible_keys = ['T'+''.join(map(str, k)) for k in product([0,1], repeat=3)]

输出:

['T000', 'T001', 'T010', 'T011', 'T100', 'T101', 'T110', 'T111']

可以结合使用字典get方法:

## as comprehension
{k: out.get(k, 0) for k in possible_keys}
## as classical loop
for key in possible_keys:
out.get(key, 0) # default value is 0

完整例子:

from itertools import product
a = [0,1,1,0]
out = {'T'+''.join(map(str, k)): v for k,v in Counter(windowed(a, n=3)).items()}
possible_keys = ['T'+''.join(map(str, k)) for k in product([0,1], repeat=3)]
out = {k: out.get(k, 0) for k in possible_keys}

输出:

{'T000': 0,
'T001': 0,
'T010': 0,
'T011': 1,
'T100': 0,
'T101': 0,
'T110': 1,
'T111': 0}

字典更新(python≥3.9)

from itertools import product
default = {'T'+''.join(map(str, k)): 0 for k in product([0,1], repeat=3)}
out = default | out # update default with the values in out