PyMC3 Beta-Binomial fails to converge on actual parameter values 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)


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)


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)



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.

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.
