我有一个二维概率矩阵。我想用二项分布从这个矩阵中抽样n次。使用numpy / scipy最优雅的方式是什么？

import numpy as np

test_mat = np.random.beta(10, 2, size=(10,3))
trials = 5

# this works but only 1 trial
np.random.binomial(1, test_mat)

# this gives error
np.random.binomial(1, test_mat, size = (10, 3, trials))

ValueError: shape mismatch: objects cannot be broadcast to a single shape
或者是我唯一的选择，首先将test_mat数组重新成形为三维（本质上，在第一和第二维的给定行和列中，为第三维的所有5列插入相同的值）。