Modeling bivariate beta distributions

Questions
3

Thank you for your answer.

It seems to me that your model considers a different probability of success for each day in both groups.

My intention was to find distributions of just two probabilities of success each of which is static for all days: the first one (X) generates successes from trials in group a, and the second one (Y) - in group b.
This approach was inspired by this series of posts and this post with pymc3 implementation.

The simplest model showing what I’m trying to achieve is the following:

with pm.Model() as model_simple:
    beta_a = pm.Beta('beta_a', alpha=PRIOR_ALPHA, beta=PRIOR_BETA)
    beta_b = pm.Beta('beta_b', alpha=PRIOR_ALPHA, beta=PRIOR_BETA)
    
    obs_a = pm.Binomial('obs_a', n=df['trials_a'].values, p=beta_a, observed=df['successes_a'].values)
    obs_b = pm.Binomial('obs_b', n=df['trials_b'].values, p=beta_b, observed=df['successes_b'].values)

    step = pm.Metropolis()
    trace_simple = pm.sample(200000, tune=20000, step=step, init='advi')

As you can see, there is only two betas, not a one for each day in each group.
The only disadvantage of this model is that betas are not correlated, which is why I’m trying to implement the model from the paper.