Causal Inference I(h) — Bayesian A/B Testing and Multi-Armed Bandits

Posteriors, expected loss, and adaptive allocation — the decision-theoretic view of experimentation

The frequentist notebooks asked "can we reject the null?" The Bayesian view reframes experimentation as a decision under uncertainty: given the data, what is the probability that B beats A, how much do we expect to lose if we ship the wrong arm, and — if we can adapt during the test — how do we send traffic to the better arm as we learn? This is how much of modern industry experimentation actually operates (VWO, Google Optimize, and every bandit-driven recommender).

We build three things:

1. Bayesian A/B — the Beta-Binomial conjugate model for conversion rates: posteriors for each arm, the probability B beats A, the expected loss of a ship decision, and how these differ from a p-value.
2. The honest caveat — Bayesian posteriors do not magically license unlimited peeking; a naive "stop when P(B>A) > 0.95" rule still over-declares winners under the null. The clean guarantee is decision-theoretic (expected loss), not a frequentist error rate.
3. Multi-armed bandits — Thompson sampling: instead of fixing a 50/50 split, allocate traffic in proportion to each arm's posterior probability of being best, minimizing regret (lost conversions from showing inferior arms).

Simulation-first with known conversion rates. Python leads (conjugate updating, expected loss, Thompson sampling from scratch); the R companion mirrors it. This links the experimentation arc to the Bayesian project set — the Beta-Binomial conjugacy here is the same machinery as the baseball shrinkage and binomial-GLMM notebooks.

1. Bayesian A/B — posteriors, P(B beats A), and expected loss

For a conversion rate, the Beta-Binomial model is conjugate: with a $\text{Beta}(1,1)$ (uniform) prior and $c$ conversions in $n$ trials, the posterior is $\text{Beta}(1+c,\,1+n-c)$. From the two arms' posteriors we read off decision-relevant quantities directly:

• $\Pr(p_B>p_A)$ — the probability B is better (Monte Carlo over posterior draws);
• the expected loss of shipping B, $\mathbb{E}[\max(p_A-p_B,0)]$ — how much conversion rate we expect to forgo if B is actually worse — and symmetrically for A. The decision rule: ship the arm whose expected loss falls below a small threshold (a "caliper" in conversion units), which bounds the downside of being wrong rather than controlling a hypothetical error rate.

On a borderline case, $\Pr(p_B>p_A)\approx0.91$ and B's expected loss is tiny — a shippable but not certain result, and the numbers say exactly how uncertain.

import numpy as np, matplotlib.pyplot as plt, warnings
from scipy.stats import beta as Beta
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
rng=np.random.default_rng(2)
pA,pB=0.10,0.112; nA=nB=1500
cA=rng.binomial(nA,pA); cB=rng.binomial(nB,pB)
aA,bA=1+cA,1+nA-cA; aB,bB=1+cB,1+nB-cB
S=200000; dA=rng.beta(aA,bA,S); dB=rng.beta(aB,bB,S)
pBA=(dB>dA).mean(); loss_B=np.mean(np.maximum(dA-dB,0)); loss_A=np.mean(np.maximum(dB-dA,0)); uplift=np.mean(dB-dA)
print(f"observed: A {cA}/{nA} = {cA/nA:.3f},  B {cB}/{nB} = {cB/nB:.3f}")
print(f"  P(p_B > p_A)            = {pBA:.3f}")
print(f"  expected uplift E[pB-pA] = {uplift:+.4f}  (95% credible [{np.percentile(dB-dA,2.5):+.4f}, {np.percentile(dB-dA,97.5):+.4f}])")
print(f"  expected loss ship B    = {loss_B:.5f}   ship A = {loss_A:.5f}")
print(f"  decision (ship B if loss < 0.0005): {'SHIP B' if loss_B<0.0005 else 'keep testing'}")
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
xg=np.linspace(0.07,0.15,400)
ax[0].plot(xg,Beta.pdf(xg,aA,bA),color=BLUE,lw=2,label=f"A posterior ({cA/nA:.3f})")
ax[0].plot(xg,Beta.pdf(xg,aB,bB),color=GREEN,lw=2,label=f"B posterior ({cB/nB:.3f})")
ax[0].fill_between(xg,Beta.pdf(xg,aA,bA),alpha=.2,color=BLUE); ax[0].fill_between(xg,Beta.pdf(xg,aB,bB),alpha=.2,color=GREEN)
ax[0].set_xlabel("conversion rate"); ax[0].set_ylabel("posterior density"); ax[0].set_title(f"Posteriors overlap -> P(B>A)={pBA:.2f}, not certainty"); ax[0].legend(fontsize=8)
ax[1].hist(dB-dA,bins=80,color=PURP,alpha=.7); ax[1].axvline(0,color=RED,lw=2)
ax[1].set_xlabel("posterior of uplift pB - pA"); ax[1].set_ylabel("draws"); ax[1].set_title(f"P(uplift>0)={pBA:.2f}; most mass positive but some below 0")
plt.tight_layout(); plt.show()
print("The Bayesian output is a full distribution of the uplift, not a yes/no verdict: it quantifies how likely B wins and how")
print("much we'd lose if we're wrong -- directly usable for a ship decision, unlike a bare p-value.")

observed: A 153/1500 = 0.102,  B 176/1500 = 0.117
  P(p_B > p_A)            = 0.911
  expected uplift E[pB-pA] = +0.0154  (95% credible [-0.0070, +0.0378])
  expected loss ship B    = 0.00047   ship A = 0.01583
  decision (ship B if loss < 0.0005): SHIP B

The Bayesian output is a full distribution of the uplift, not a yes/no verdict: it quantifies how likely B wins and how
much we'd lose if we're wrong -- directly usable for a ship decision, unlike a bare p-value.

2. The honest caveat — Bayesian posteriors don't license free peeking

A common claim is that Bayesian A/B "solves" the peeking problem. It does not, in the frequentist sense. If you monitor continuously and stop the moment $\Pr(p_B>p_A)>0.95$, then even when the two arms are identical you will cross that threshold by chance surprisingly often — the posterior wanders just as the running z-statistic did (notebook 1d). We simulate an A/A test under this rule and find it declares a "winner" far more than 5% of the time.

What Bayesian decision theory does give you is a clean decision-theoretic guarantee: stopping when the expected loss of shipping the leader is below a small caliper bounds how much you can expect to lose, regardless of when you stop. That is a statement about decision quality, not about a long-run false-positive rate — a different (and often more useful) contract, but one worth stating honestly rather than overselling.

def bayes_peek(seed, pA=0.10, pB=0.10, N=4000, thr=0.95, step=50, draws=3000):
    r=np.random.default_rng(seed); a=np.cumsum(r.random(N)<pA); b=np.cumsum(r.random(N)<pB)
    for n in range(step,N+1,step):
        d1=r.beta(1+a[n-1],1+n-a[n-1],draws); d2=r.beta(1+b[n-1],1+n-b[n-1],draws)
        pp=(d2>d1).mean()
        if pp>thr or pp<1-thr: return True
    return False
false_win=np.mean([bayes_peek(s) for s in range(600)])
fixed_win=np.mean([ (lambda s: (lambda r,a,b: (r.beta(1+a,1+4000-a,3000)> r.beta(1+b,1+4000-b,3000)).mean())(np.random.default_rng(s), np.random.default_rng(s).binomial(4000,0.10), np.random.default_rng(1000+s).binomial(4000,0.10)))(s)>0.95 for s in range(600)])
print(f"A/A test (arms identical), 'stop when P(B>A)>0.95' rule:")
print(f"  continuous peeking declares a winner {false_win:.2f} of the time  <- NOT 5%: peeking still inflates error")
print(f"  (a single fixed-n look at P>0.95 fires ~{fixed_win:.2f} -- as designed)")
fig,ax=plt.subplots(figsize=(7.5,4))
ax.bar(["fixed single look","continuous peeking\n(stop at P>0.95)"],[fixed_win,false_win],color=[GREEN,RED])
ax.axhline(0.05,color="k",ls="--",label="5% reference"); ax.set_ylabel("P(declare a winner) under the null"); ax.set_title("Bayesian posteriors do not make peeking free"); ax.legend()
for i,v in enumerate([fixed_win,false_win]): ax.text(i,v+0.01,f"{v:.2f}",ha="center")
plt.tight_layout(); plt.show()
print("Use Bayesian quantities for the DECISION (expected loss below a caliper), not as a license to stop the instant a")
print("probability threshold is crossed -- that reintroduces the peeking inflation of notebook 1d.")

A/A test (arms identical), 'stop when P(B>A)>0.95' rule:
  continuous peeking declares a winner 0.58 of the time  <- NOT 5%: peeking still inflates error
  (a single fixed-n look at P>0.95 fires ~0.06 -- as designed)

Use Bayesian quantities for the DECISION (expected loss below a caliper), not as a license to stop the instant a
probability threshold is crossed -- that reintroduces the peeking inflation of notebook 1d.

3. Multi-armed bandits — Thompson sampling and regret

A/B testing splits traffic fixed at 50/50 (or 1/K) for the whole test, so half your users keep seeing the worse arm until it ends — that lost conversion is regret. A multi-armed bandit adapts: it sends more traffic to arms that look better as data accrues, trading exploration for exploitation. Thompson sampling is the elegant Bayesian rule — at each step, draw one sample from each arm's posterior and play the arm with the highest draw; arms that are probably-best get most traffic, but uncertain arms still get explored in proportion to their posterior chance of winning.

Across four arms with true rates 0.10–0.13, we compare cumulative regret of uniform A/B, epsilon-greedy, and Thompson sampling. Thompson achieves the lowest regret — the same experiment, far fewer conversions sacrificed to inferior arms. The trade-off: bandits optimize earnings during the test rather than clean end-of-test inference, so they suit ongoing optimization (which of these headlines?) more than a one-off causal readout.

true=np.array([0.10,0.11,0.12,0.13]); K=len(true); T=20000; best=true.max()
def run_AB(seed):
    r=np.random.default_rng(seed); hist=np.empty(T)
    for t in range(T): a=t%K; hist[t]=true[a]; _=r.random()<true[a]
    return hist
def run_eps(seed,eps=0.1):
    r=np.random.default_rng(seed); pulls=np.ones(K); rew=np.full(K,0.5); hist=np.empty(T)
    for t in range(T):
        a=r.integers(K) if r.random()<eps else int(np.argmax(rew/pulls))
        x=r.random()<true[a]; pulls[a]+=1; rew[a]+=x; hist[t]=true[a]
    return hist
def run_ts(seed):
    r=np.random.default_rng(seed); a_=np.ones(K); b_=np.ones(K); hist=np.empty(T); alloc=np.zeros(K)
    for t in range(T):
        a=int(np.argmax(r.beta(a_,b_))); x=r.random()<true[a]; a_[a]+=x; b_[a]+=1-x; hist[t]=true[a]; alloc[a]+=1
    return hist, alloc
regret=lambda h: np.cumsum(best-h)
R_AB=np.mean([regret(run_AB(s))[-1] for s in range(60)])
R_eps=np.mean([regret(run_eps(s))[-1] for s in range(60)])
tss=[run_ts(s) for s in range(60)]; R_ts=np.mean([regret(h)[-1] for h,_ in tss])
alloc=np.mean([a for _,a in tss],axis=0)
print(f"Cumulative regret over T={T:,} pulls (lost conversions vs always-best, lower is better):")
print(f"  uniform A/B     : {R_AB:7.1f}")
print(f"  epsilon-greedy  : {R_eps:7.1f}")
print(f"  Thompson sampling: {R_ts:6.1f}   <- lowest regret")
print(f"  Thompson traffic allocation by arm (rates {true}): {np.round(alloc/T,2)}  -> most traffic to the best arm")
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].plot(regret(run_AB(0)),color=GREY,lw=2,label="uniform A/B")
ax[0].plot(regret(run_eps(0)),color=ORANGE,lw=2,label="epsilon-greedy")
ax[0].plot(regret(run_ts(0)[0]),color=GREEN,lw=2,label="Thompson")
ax[0].set_xlabel("users"); ax[0].set_ylabel("cumulative regret"); ax[0].set_title("Thompson sampling accrues the least regret"); ax[0].legend()
ax[1].bar(range(K),alloc/T,color=[GREY,GREY,GREY,GREEN]); ax[1].set_xticks(range(K)); ax[1].set_xticklabels([f"{r:.2f}" for r in true])
ax[1].set_xlabel("arm true rate"); ax[1].set_ylabel("share of traffic"); ax[1].set_title("Thompson concentrates traffic on the best arm")
plt.tight_layout(); plt.show()
print("Thompson sends most traffic to the best arm while still exploring -- minimizing lost conversions. Bandits optimize")
print("cumulative reward during the test, the right objective for continuous optimization (vs a clean one-shot causal estimate).")

Cumulative regret over T=20,000 pulls (lost conversions vs always-best, lower is better):
  uniform A/B     :   300.0
  epsilon-greedy  :   102.2
  Thompson sampling:   94.6   <- lowest regret
  Thompson traffic allocation by arm (rates [0.1  0.11 0.12 0.13]): [0.03 0.09 0.21 0.68]  -> most traffic to the best arm

Thompson sends most traffic to the best arm while still exploring -- minimizing lost conversions. Bandits optimize
cumulative reward during the test, the right objective for continuous optimization (vs a clean one-shot causal estimate).

4. Summary

The Bayesian view turns experimentation into a decision under uncertainty:

• Bayesian A/B (Beta-Binomial conjugate) reports the full posterior of the uplift — $\Pr(p_B>p_A)$, a credible interval, and the expected loss of a ship decision — which is directly actionable in a way a p-value is not; on the borderline case, $\Pr(B>A)=0.91$ with a small expected loss said "probably ship, here's the risk."
• The honest caveat: Bayesian posteriors do not make peeking free — a naive rule that stops the moment P(B beats A) exceeds 0.95 over-declared winners under the null. The Bayesian guarantee is decision-theoretic (bounded expected loss), not a frequentist error rate; for strict error control, use the always-valid methods of notebook 1d.
• Multi-armed bandits / Thompson sampling adaptively route traffic to the better arm, cutting regret roughly threefold versus a uniform A/B — the right tool when the goal is to earn while learning rather than produce a single clean causal estimate.

Guidance: use Bayesian A/B for interpretable, decision-focused readouts (uplift distribution + expected loss with a caliper); use bandits for ongoing optimization across many variants; and reach for frequentist always-valid methods (1d) when you need a guaranteed false-positive rate. Cross-links: the Beta-Binomial conjugacy is the exact machinery of the Bayesian arc's baseball shrinkage and binomial-GLMM notebooks; Thompson sampling is Bayesian decision theory in action (posterior draws → action), connecting to the priors/posteriors thread throughout that arc; and the peeking caveat ties back to Sequential Testing (1d). The R companion reproduces the conjugate updating, expected loss, and Thompson sampling in base R.