Bayesian Models and Their Role in A/B Testing (R Notes)
Two versions of a landing page are live at once — Page A and Page B — and the question the marketing team keeps asking is simple: which one converts better, and how sure are we? A frequentist test gives a yes/no answer once a fixed sample size is hit. A Bayesian model gives something more useful along the way: a full probability distribution over "what the true conversion rate might be," updated continuously as new visitors arrive, and a direct answer to "what's the probability B is actually better than A?"
• A Bayesian model starts with a prior belief about the unknown quantity.
• It updates that belief using observed data through Bayes' theorem.
• The result is a posterior distribution — not a single number, but a full range of plausible values with their relative likelihoods attached.
Core Idea: Bayes' Theorem
Everything here comes from one equation:
In plain terms: start with what you believed before seeing any data (the prior), weigh it by how well each possible value explains the data you actually observed (the likelihood), and the result — after rescaling so it's a proper probability distribution — is your updated belief (the posterior).
For a conversion rate problem, there's a convenient shortcut called a conjugate prior. If the prior for a conversion rate is a Beta distribution, and the data is a count of successes and failures (a Binomial process), the posterior is also a Beta distribution, with parameters updated by simple addition:
The Dataset: 10 Stores, One Marketing Campaign
Instead of just two pages, here's a more realistic business analytics scenario: a company ran the same promotional campaign across 10 stores and logged visitor and conversion counts for each. The question isn't just "which store did best on paper" — it's "which store is most likely to actually be the best performer, once sampling noise is accounted for?" A store with 150 visitors and a 16.7% conversion rate is a more trustworthy signal than one with 20 visitors and a lucky 20% rate.
| Store_ID | Region | Visitors | Conversions | Conversion_Rate |
|---|---|---|---|---|
| S01 | North | 120 | 14 | 11.7% |
| S02 | South | 95 | 9 | 9.5% |
| S03 | East | 140 | 22 | 15.7% |
| S04 | West | 110 | 11 | 10.0% |
| S05 | North | 130 | 19 | 14.6% |
| S06 | South | 100 | 8 | 8.0% |
| S07 | East | 150 | 25 | 16.7% |
| S08 | West | 90 | 7 | 7.8% |
| S09 | North | 125 | 17 | 13.6% |
| S10 | South | 105 | 10 | 9.5% |
Five columns: Store_ID, Region, Visitors, Conversions, and the derived Conversion_Rate. Ranking by that last column alone is exactly the trap a Bayesian model avoids — S03's 15.7% and S07's 16.7% look close, but they're built on different sample sizes, and a raw percentage doesn't say how much to trust either number.
Worked Example: Comparing Two Landing Pages
Before scaling up to all 10 stores, the same Beta-Binomial mechanics from earlier still apply — just to a smaller pair. Both start with no prior opinion about which converts better — a flat, uninformative prior, Beta(1, 1), which just says every conversion rate between 0 and 1 is equally plausible before any data comes in.
| Page | Visitors | Conversions | Prior | Posterior |
|---|---|---|---|---|
| A | 100 | 8 | Beta(1, 1) | Beta(1+8, 1+92) = Beta(9, 93) |
| B | 100 | 12 | Beta(1, 1) | Beta(1+12, 1+88) = Beta(13, 89) |
What the posteriors say
Page B's posterior mean sits higher, but the two credible intervals overlap — so the raw numbers alone don't settle it. The Bayesian model's real advantage is answering the actual question directly: what's the probability that B's true conversion rate is higher than A's? That's done by drawing a large number of random samples from each posterior and simply counting how often B's sample beats A's.
Scaling Up: Ranking All 10 Stores
The same update rule applies to every row of the dataset independently — each store gets its own posterior. With 10 stores instead of 2, "which one is better" becomes "which one is most likely the best out of all ten," which needs one extra step: draw a sample from all 10 posteriors at once, note which store wins that particular draw, and repeat that thousands of times. The fraction of draws each store wins is its probability of being the best performer.
| Store_ID | Region | Posterior mean | 95% credible interval | P(best performer) |
|---|---|---|---|---|
| S01 | North | 12.3% | [7.1%, 18.7%] | 3.7% |
| S02 | South | 10.3% | [5.1%, 17.0%] | 1.4% |
| S03 | East | 16.2% | [10.6%, 22.7%] | 26.3% |
| S04 | West | 10.7% | [5.7%, 17.0%] | 1.4% |
| S05 | North | 15.2% | [9.6%, 21.7%] | 16.8% |
| S06 | South | 8.8% | [4.2%, 15.0%] | 0.4% |
| S07 | East | 17.1% | [11.6%, 23.5%] | 38.0% |
| S08 | West | 8.7% | [3.9%, 15.2%] | 0.5% |
| S09 | North | 14.2% | [8.7%, 20.7%] | 10.5% |
| S10 | South | 10.3% | [5.3%, 16.7%] | 1.1% |
S07 comes out ahead with a 38% probability of being the true best performer — notably higher than its raw 16.7% conversion rate alone would suggest, because that rate is backed by the largest sample in the dataset (150 visitors). S03 is a close second at 26%, despite a slightly lower raw conversion rate, since its sample size (140 visitors) also earns it a tighter, more trustworthy posterior. Compare that to S06 and S08, whose posterior means look unremarkable and whose probability of being the best sits below 1% — small samples with middling rates just don't have a realistic shot at the top spot once uncertainty is priced in.
The R Code, Explained Line by Line
Here is the full script — extended to loop over all 10 stores in the dataset — followed by a breakdown of what each part is doing and why.
# --- Bayesian ranking of 10 stores: Beta-Binomial model ---
set.seed(42)
# The business analytics dataset: 10 records, 5 columns
stores <- data.frame(
Store_ID = c("S01","S02","S03","S04","S05","S06","S07","S08","S09","S10"),
Region = c("North","South","East","West","North",
"South","East","West","North","South"),
Visitors = c(120, 95, 140, 110, 130, 100, 150, 90, 125, 105),
Conversions = c(14, 9, 22, 11, 19, 8, 25, 7, 17, 10)
)
stores$Conversion_Rate <- round(stores$Conversions / stores$Visitors * 100, 1)
# Uninformative prior: Beta(1, 1) is uniform over [0, 1]
prior_alpha <- 1
prior_beta <- 1
# Posterior parameters for every store at once (vectorized, no loop needed)
stores$post_alpha <- prior_alpha + stores$Conversions
stores$post_beta <- prior_beta + (stores$Visitors - stores$Conversions)
stores$post_mean <- stores$post_alpha / (stores$post_alpha + stores$post_beta)
# 95% credible interval per store
ci <- mapply(function(a, b) qbeta(c(0.025, 0.975), a, b),
stores$post_alpha, stores$post_beta)
stores$ci_low <- round(ci[1, ], 4)
stores$ci_high <- round(ci[2, ], 4)
print(stores[, c("Store_ID","Region","post_mean","ci_low","ci_high")])
# Monte Carlo: draw from all 10 posteriors at once, find the winner per draw
n_samples <- 100000
sample_matrix <- sapply(1:nrow(stores), function(i)
rbeta(n_samples, stores$post_alpha[i], stores$post_beta[i]))
colnames(sample_matrix) <- stores$Store_ID
winner_per_draw <- apply(sample_matrix, 1, which.max)
prob_best <- table(factor(winner_per_draw, levels = 1:nrow(stores))) / n_samples
names(prob_best) <- stores$Store_ID
cat("Probability each store is the best performer:\n")
print(round(prob_best, 4))
cat("\nTop store:", names(which.max(prob_best)), "\n")
# Visualize all 10 posteriors together
plot(NULL, xlim = c(0, 0.3), ylim = c(0, 25),
xlab = "conversion rate", ylab = "density",
main = "Posterior distributions: all 10 stores")
colors <- rainbow(nrow(stores))
for (i in 1:nrow(stores)) {
lines(density(sample_matrix[, i]), col = colors[i], lwd = 2)
}
legend("topright", legend = stores$Store_ID, col = colors, lwd = 2, cex = 0.7)
1 S01 North 0.1230 0.0711 0.1865
2 S02 South 0.1031 0.0511 0.1704
3 S03 East 0.1620 0.1063 0.2267
4 S04 West 0.1071 0.0571 0.1704
5 S05 North 0.1515 0.0958 0.2173
6 S06 South 0.0882 0.0416 0.1502
7 S07 East 0.1711 0.1156 0.2345
8 S08 West 0.0870 0.0387 0.1521
9 S09 North 0.1417 0.0869 0.2072
10 S10 South 0.1028 0.0530 0.1666
Probability each store is the best performer:
S01 S02 S03 S04 S05 S06 S07 S08 S09 S10
0.0367 0.0143 0.2625 0.0142 0.1676 0.0039 0.3801 0.0047 0.1048 0.0112
Top store: S07
set.seed(42) — fixes R's random number generator so the Monte Carlo simulation later in the script gives the same result every time it's re-run.data.frame(Store_ID = ..., Region = ..., Visitors = ..., Conversions = ...) — loads the 10-record, 5-column dataset directly into R's native tabular structure, the same way it would come in from a CSV in a real analytics pipeline.stores$Conversion_Rate <- round(...) — derives the fifth column from the other four. This is the number a plain spreadsheet ranking would use — and exactly the number the rest of the script shows isn't trustworthy on its own.prior_alpha <- 1; prior_beta <- 1 — the same flat, uninformative starting belief used for the two-page example, applied identically to every store.stores$post_alpha <- prior_alpha + stores$Conversions — this line updates all 10 stores at once. Because R operates on whole columns by default, there's no loop needed to apply the Beta-Binomial update rule to every row.mapply(function(a, b) qbeta(...), ...) — runs the credible-interval calculation once per store. mapply is R's way of applying a function across two parallel vectors (here, each store's α and β) instead of writing a manual for-loop.sapply(1:nrow(stores), function(i) rbeta(...)) — draws 100,000 samples from each of the 10 posteriors, arranging the result as a matrix with one column per store. This is the multi-store generalization of the two-sample rbeta() call from the earlier example.apply(sample_matrix, 1, which.max) — for every one of the 100,000 simulated draws, finds which store had the highest sampled conversion rate in that draw. Running this row by row across the matrix is what turns 10 separate distributions into a single ranked competition.table(factor(winner_per_draw, ...)) / n_samples — counts how many of the 100,000 draws each store won, then converts that count into a probability. This is the direct answer to "which store is really the best," accounting for how much (or little) data backs each one.Reading the output: S07 wins the head-to-head simulation most often, backing up the 38% "probability of being best" figure from the ranking table earlier. That's a meaningfully strong signal — not certain enough to shut down every other store, but strong enough that most teams would route more budget toward S07 and S03 while continuing to collect data from the rest.
Why R Specifically
Bayesian modeling can be done in Python, Julia, or by hand — but R has a few properties that make it a natural fit for this kind of work:
- → Every common distribution ships built in.
dbeta,pbeta,qbeta, andrbeta(density, cumulative probability, quantile, and random draws) exist for the Beta distribution and dozens of others, with no extra installation. - → Vectorized by default.
mean(samples_B > samples_A)compares 100,000 pairs of numbers and reduces them to one probability in a single line, with no explicit loop. - → Built for statistical output first. Base R plotting (
density(),plot()) and summary functions were designed around exactly this kind of "show me the distribution, not just a point estimate" workflow. - → A mature Bayesian ecosystem sits on top of base R — packages like
rstanarm,brms, andbayesmextend this same conjugate-model idea to full hierarchical and regression-based Bayesian models when a simple closed-form update like this one isn't enough.
Summary Table
| Aspect | Detail |
|---|---|
| Purpose | Estimate an unknown quantity (here, a store's conversion rate) as a full probability distribution, not a single number |
| Dataset | 10 records × 5 columns: Store_ID, Region, Visitors, Conversions, Conversion_Rate |
| Prior | Beta(1, 1) — uninformative, flat over [0, 1] |
| Update rule | Beta(α, β) + k successes, n trials → Beta(α+k, β+n−k) |
| Posterior mean | Closed form: α / (α + β) |
| Credible interval | qbeta() — exact, no simulation required |
| Ranking many posteriors | Monte Carlo: sample all, count how often each one wins per draw |
| R functions used | qbeta, rbeta, sapply, mapply, apply, table, density, plot |
Advantages
- → Produces a direct probability statement ("82% chance B is better") instead of a binary reject/fail-to-reject verdict.
- → Can be updated continuously as new data arrives — no need to wait for a fixed sample size before looking at results.
- → Naturally expresses uncertainty through the width of the posterior, rather than needing a separate confidence-interval calculation bolted on afterward.
Disadvantages
- → Choice of prior can influence results, especially with small sample sizes — a badly chosen prior can bias early conclusions.
- → Closed-form posteriors like Beta-Binomial only exist for certain prior/likelihood pairings; more complex models often require simulation-heavy methods like MCMC.
- → Continuously checking results ("peeking") still needs care — a high probability early on can be a small-sample fluke that a few more visitors quietly correct.