Product of two probabilities distributions

[Marion2025]

You might be able to just write a custom distribution. These two posts might be helpful.

Sum of Product of Bernoulli and Normal Random Variables

The product of a normal and Rademacher variables, independent from each other

There also seems to be a similar post on this where a custom distribution and potential were used

Likelihood given by product of two distributions?

I currently have a simple model: model = pm.Model() with model: # Priors for unknown model parameters hyper_mu = pm.Normal('hyper_mu', mu=0, sd=10) hyper_sd = pm.Gamma('hyper_sd', alpha=0.01, beta=0.001) c = pm.Normal('c', mu=hyper_mu, sd=hyper_sd, shape=(4, 4)) # Expected probabilities (4 x 4) p = sigmoid(c) # Likelihood (sampling distribution) of observations. data, N_data = (4 x 4) pm.Binomial('L', n=N_data, p=p, observed=data_successes) Where I want to obtain 4x4 param…