Modeling bivariate beta distributions

Questions
6

Okay, maybe I should have started with a simpler model instead of a bivariate beta. Beta can be approximated with normal so let’s replace multivariate beta with multivariate normal for simplicity:

prior_mean = PRIOR_ALPHA / (PRIOR_ALPHA + PRIOR_BETA)
with pm.Model() as normal_model:
    prob_of_success = pm.MvNormal(
        'prob_of_success',
        mu=np.array([prior_mean, prior_mean]), 
        cov=np.array([[0.002, 0.0005], [0.0005, 0.002]]),
        shape=2
    )
    
    obs = pm.Binomial(
        'obs',
        n=df[['trials_a', 'trials_b']].values,
        p=prob_of_success,
        observed=df[['successes_a', 'successes_b']].values
    )

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

This model shows no correlation of probabilities of success, too. Why?

np.corrcoef(trace_normal.get_values('prob_of_success').T)
# array([[1.        , 0.00471091],
#        [0.00471091, 1.        ]])

I am really missing something…