Efficient way to sample from matrix of probabilities n times with binomial distribution

I have a 2-dim matrix of probabilities. I want to sample using a binomial distribution from this matrix n times. What is the most elegant way to do this using 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

Or is my only option to first reshape the test_mat array to be 3 dimensional (in essence slotting in the same values for all 5 columns in the 3rd dimension from a given row and column in the 1st and 2nd dims).