bayesian – Troubleshooting Convergence in PyMC3 Beta-Binomial Models

bayesianbeta distributionbinomial distributionpymc

Something is not performing as expected with PyMC. I'm trying a simple Beta-Binomial conjugate prior model, trying to recover known parameters.

Control data

from scipy.stats import beta
import numpy as np
import seaborn as sns
def find_ab(x, search=range(2, 22)):
    avg, score, best_a, best_b = 1000, 1000, 1, 1
    for i in search:
        for j in search:
            samples = beta.rvs(i, j, size=1000)
            new_avg = np.mean(samples)
            new_score = abs(new_avg - x)
            if new_score < score:
                best_a, best_b = i,j
                avg = new_avg
                score = new_score
    return best_a, best_b, avg

ctrl = 0.61
a_ctrl, b_ctrl ,avg_ctrl = find_ab(ctrl)
print(f"a={a_ctrl}, b={b_ctrl}, avg={avg_ctrl}")
samples_ctrl = beta.rvs(a_ctrl, b_ctrl, size=10_000)
sns.kdeplot(samples_ctrl)

ctrl

Test data

test = ctrl + 0.12
a_test, b_test ,avg_test = find_ab(test)
print(f"a={a_test}, b={b_test}, avg={avg_test}")
samples_test = beta.rvs(a_test, b_test, size=10_000)
sns.kdeplot(samples_test)

test data

As you can see, control is centered on 0.61 and test on 0.73. I'll generate samples using Numpy, effectively 10,000 independent Bernoulli trials, each, with 0.61 and 0.73 as inputs, respectively for p.

ris_ctrl = np.random.binomial(n=1,p=avg_ctrl,size=10_000)
ris_test = np.random.binomial(n=1,p=avg_test,size=8_000)

And the modeling in PyMC3

import pymc as pm
with pm.Model() as model_ctrl:
    theta_ctrl = pm.Beta('theta', a_ctrl, b_ctrl)
    lik_ctrl = pm.Binomial('likelihood',p=theta_ctrl, n=len(ris_ctrl), observed=ris_ctrl)
    trace_ctrl = pm.sample(chains=4, draws=4000)
posterior_ctrl = trace_ctrl.posterior['theta'].stack(sample=("chain", "draw")).values
pm.plot_trace(trace_ctrl, compact=False)

with pm.Model() as model_test:
    theta_test = pm.Beta('theta', a_test, b_test)
    lik_test = pm.Binomial('likelihood',p=theta_test, n=len(ris_test), observed=ris_test)
    trace_test = pm.sample(chains=4, draws=4000)
posterior_test = trace_test.posterior['theta'].stack(sample=("chain", "draw")).values
pm.plot_trace(trace_test, compact=False)

ctrl post

ctrl test

But when I look at the traces plots, the ctrl and test have centered on 0.00061 and 0.00091, respectively. I have no idea why the means are magnitudes lower than the parameters that they should recover.

What is going on here? Is this an issue with how I've defined prior and likelihood or an issue with PyMC3?

I know that Beta-Binomial is a conjugate prior so this is unnecessary; but this is effectively why I chose this prior-likelihood as a test use case.

Best Answer

You specify the number of trials parameter, $n$, inconsistently.

First you simulate Bernoulli random variables with np.random.binomial(n=1, p=p, size=size). Then in the likelihood you have pm.Binomial("likelihood", p=theta, n=len(sample), observed=sample).

Replace n=len(sample) with n=1 to get the expected result.

You can do a back-of-the-envelope calculation to understand why you got the estimates that you report: 0.00061 ≈ 0.61 / 10_000 and 0.00091 ≈ 0.73 / 8_000.

Related Question