Thanks for the quick reply as always @junpenglao.
Unfortunately it did not work…
I have used L = pm.Bernoulli('L', p=p, observed=bernoulli_data_parameter)
where bernoulli_data_parameter.shape = (4, 4, 1000).
The error:
ValueError: Input dimension mis-match. (input[0].shape[2] = 1000, input[1].shape[2] = 4)