根据先前的值对numpy代码进行矢量化处理 [英] Vectorize numpy code with operation depending on previous value
问题描述
以下代码对一个可以在任何时间采样3个不同状态的系统进行建模,并且这些状态之间的恒定转换概率由矩阵prob_nor
给出.因此,trace
中的每个点都取决于先前的状态.
The following code models a system that can sample 3 different states at any time, and the constant transition probability between those states is given by the matrix prob_nor
. Threrefore, each point in trace
depends on the previous state.
n_states, n_frames = 3, 1000
state_val = np.linspace(0, 1, n_states)
prob = np.random.randint(1, 10, size=(n_states,)*2)
prob[np.diag_indices(n_states)] += 50
prob_nor = prob/prob.sum(1)[:,None] # transition probability matrix,
# row sum normalized to 1.0
state_idx = range(n_states) # states is a list of integers 0, 1, 2...
current_state = np.random.choice(state_idx)
trace = []
sigma = 0.1
for _ in range(n_frames):
trace.append(np.random.normal(loc=state_val[current_state], scale=sigma))
current_state = np.random.choice(state_idx, p=prob_nor[current_state, :])
以上代码中的循环使其运行非常慢,尤其是当我必须对数百万个数据点进行建模时.有什么方法可以向量化/加速它?
The loop in the above code makes it run pretty slow, specially when I have to model millions of data points. Is there any way to vectorize/accelerate it?
推荐答案
尽快卸载概率计算:
possible_paths = np.vstack(
np.random.choice(state_idx, p=prob_nor[curr_state, :], size=n_frames)
for curr_state in range(n_states)
)
然后,您可以简单地按照路径进行查找:
Then you can simply do a lookup to follow your path:
path_trace = [None]*n_frames
for step in range(n_frames):
path_trace[step] = possible_paths[current_state, step]
current_state = possible_paths[current_state, step]
一旦有了路径,就可以计算出轨迹:
Once you have your path, you can compute your trace:
sigma = 0.1
trace = np.random.normal(loc=state_val[path_trace], scale=sigma, size=n_frames)
比较时间:
纯python for
循环
Pure python for
loop
%%timeit
trace_list = []
current_state = np.random.choice(state_idx)
for _ in range(n_frames):
trace_list.append(np.random.normal(loc=state_val[current_state], scale=sigma))
current_state = np.random.choice(state_idx, p=prob_nor[current_state, :])
结果:
30.1 ms ± 436 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
矢量化查找:
%%timeit
current_state = np.random.choice(state_idx)
path_trace = [None]*n_frames
possible_paths = np.vstack(
np.random.choice(state_idx, p=prob_nor[curr_state, :], size=n_frames)
for curr_state in range(n_states)
)
for step in range(n_frames):
path_trace[step] = possible_paths[current_state, step]
current_state = possible_paths[current_state, step]
trace = np.random.normal(loc=state_val[path_trace], scale=sigma, size=n_frames)
结果:
641 µs ± 6.03 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
加速约50倍.
这篇关于根据先前的值对numpy代码进行矢量化处理的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!