Efficient Computation Across Many Treatment Arms

Why Vectorize the Comparison?

With two treatment arms, it is easy to compare posterior draws one pair at a time. With many arms, that approach becomes cumbersome: every draw must be compared across all arms, and the calculation must be repeated many times to estimate the probability that each arm is best. A matrix representation makes the repeated comparison explicit, reduces opportunities for indexing errors, and allows the same calculation to be applied to any number of arms.

The target is not simply the arm with the largest observed response rate or posterior mean. Under Thompson sampling, we draw a response rate from each arm’s posterior distribution and select the arm with the largest draw. Repeating this operation estimates how often each arm would be selected under posterior sampling, conditional on the trial data and model.

Definition: Let $K$ be the number of treatment arms, and let $\theta_k$ denote the response rate for arm $k$, for $k=1,\ldots,K$. Given the observed trial data $D$, the posterior best-arm probability for arm $k$ is $$ q_k=P\left(\theta_k=\max_{1\leq j\leq K}\theta_j\mid D\right). $$ When posterior distributions are continuous, ties have probability zero. The values $q_1,\ldots,q_K$ are nonnegative and sum to one. They describe posterior uncertainty about which arm has the largest response rate, not the predictive response probability for a future participant.

Draws as a Matrix

Suppose the posterior for the $K$ response rates factors into independent arm-specific posteriors given $D$. For each Monte Carlo replicate $m=1,\ldots,M$, draw one value from each arm’s posterior. Store the result in an $M$-by-$K$ matrix $Z$, where $Z_{mk}$ is the draw for arm $k$ in replicate $m$. Each row is one joint posterior draw across the arms; each column contains $M$ draws for one arm.

The calculation has two operations. First, find the largest entry in each row. Second, record which column contains that row maximum. Define $I_{mk}=1$ if arm $k$ has the largest draw in row $m$, and $I_{mk}=0$ otherwise. The Monte Carlo estimate of the posterior best-arm probability is the fraction of rows won by arm $k$:

\widehat{q}_k=\frac{1}{M}\sum_{m=1}^{M} I_{mk}

In vectorized terms, the row maxima are obtained by reducing across the $K$ columns, and the winning-arm indicators are computed for all rows and columns together. This avoids writing a separate comparison procedure for each arm pair or each draw. The matrix is a computational device: the statistical target remains the posterior probability $q_k$.

Assumptions for this calculation: The arm-specific posteriors are correctly specified for the data and factor independently conditional on $D$; a response rate has the same interpretation and direction across arms; and the Monte Carlo replicates are independent draws from the joint posterior. For continuous posteriors, row maxima are unique with probability one. If a model permits ties, the computation needs an explicitly defined tie rule.

Why the Empirical Winner Fraction Works

For a fixed arm $k$, each $I_{mk}$ is an indicator: it equals one when arm $k$ is best in replicate $m$, and zero otherwise. Under the stated posterior-sampling assumptions, $I_{mk}$ is a Bernoulli random variable with success probability $q_k$. Independent Monte Carlo replicates therefore give

E(\widehat{q}_k\mid D)=q_k,\qquad \operatorname{Var}(\widehat{q}_k\mid D)=\frac{q_k(1-q_k)}{M}.

The expectation follows by averaging $M$ indicators, each with expectation $q_k$. Independence gives the variance: the variance of their sum is $M q_k(1-q_k)$, and dividing by $M^2$ yields $q_k(1-q_k)/M$. Since $q_k$ is unknown, a customary estimated Monte Carlo standard error is

\widehat{SE}_{MC}(\widehat{q}_k) =\sqrt{\frac{\widehat{q}_k(1-\widehat{q}_k)}{M}}.

This standard error describes simulation error in estimating a posterior probability. It is not a measure of uncertainty in the underlying treatment response rate; that uncertainty is represented by the posterior distributions themselves.

Worked Example: Comparing Posterior Response Rates for Three Arms

Consider a hypothetical trial with binary responses, where a response is favorable and a larger response probability is better. Suppose each arm began with a $\operatorname{Beta}(4,10)$ prior. The observed data are 8 responses among 30 participants on arm A, 11 among 30 on arm B, and 9 among 30 on arm C. With independent Beta priors and Bernoulli outcomes, the arm posteriors are independent Beta distributions.

For a Beta-Binomial update, add the number of responses to the first parameter and the number of nonresponses to the second. Thus the posteriors are $\operatorname{Beta}(12,32)$ for A, $\operatorname{Beta}(15,29)$ for B, and $\operatorname{Beta}(13,31)$ for C. For example, arm A has 30 − 8 = 22 nonresponses, giving $\operatorname{Beta}(4+8,10+22)=\operatorname{Beta}(12,32)$.

For illustration, suppose eight Monte Carlo replicates from these posteriors give the following response-rate draws. These are hypothetical draws displayed to show the calculation; eight draws are far too few for a precise practical estimate.

Replicate	Arm A	Arm B	Arm C	Winner
1	0.31	0.35	0.28	B
2	0.27	0.24	0.29	C
3	0.33	0.30	0.25	A
4	0.22	0.31	0.27	B
5	0.29	0.28	0.32	C
6	0.26	0.34	0.30	B
7	0.35	0.29	0.31	A
8	0.30	0.27	0.26	A

The row maxima occur for A three times, B three times, and C twice. Dividing each count by $M=8$ gives

\widehat{q}_A=\frac{3}{8}=0.375,\qquad \widehat{q}_B=\frac{3}{8}=0.375,\qquad \widehat{q}_C=\frac{2}{8}=0.250.

These estimates sum to one. They do not establish that A and B have exactly equal best-arm probabilities: the apparent equality is based on only eight replicates. The calculation estimates how often each arm’s posterior draw is the largest, not which arm has the largest posterior mean.

Many Arms and a Gamma-Poisson Example

The matrix calculation does not depend on a Beta model. For example, suppose each arm has a Poisson event rate and larger rates are clinically desirable. If arm $k$ has total count $Y_k$ over exposure $E_k$, and its prior is Gamma with shape $a_k$ and rate $b_k$, then its posterior is Gamma with shape $a_k+Y_k$ and rate $b_k+E_k$. Here the Gamma distribution is parameterized by shape and rate, so its mean is shape divided by rate. With independent priors and arm-specific Poisson likelihoods, the posteriors factor across arms.

Worked Example: Vectorized Winners for Four Poisson-Rate Arms

In a hypothetical example, four arms each contribute 20 person-years of exposure and have observed event totals 9, 11, 10, and 8 for arms A through D. Suppose every arm has an independent $\operatorname{Gamma}(2,4)$ prior, using the shape-rate parameterization. The posterior shapes are $2+9=11$, $2+11=13$, $2+10=12$, and $2+8=10$. The posterior rate for each arm is $4+20=24$. Thus the posterior means, in events per person-year, are $11/24=0.4583$, $13/24=0.5417$, $12/24=0.5000$, and $10/24=0.4167$.

Suppose six illustrative Monte Carlo replicates produce the following rates, in events per person-year. The numbers are hypothetical posterior draws, not additional observed trial data.

Replicate	Arm A	Arm B	Arm C	Arm D	Winner
1	0.42	0.38	0.51	0.47	C
2	0.36	0.44	0.41	0.39	B
3	0.55	0.49	0.46	0.52	A
4	0.40	0.45	0.48	0.50	D
5	0.47	0.43	0.52	0.46	C
6	0.51	0.54	0.49	0.50	B

Counting the row maxima gives one win for A, two for B, two for C, and one for D. Therefore the estimates are $1/6=0.1667$, $2/6=0.3333$, $2/6=0.3333$, and $1/6=0.1667$, respectively. These four estimates sum to one. Their interpretation depends on the stated clinical direction: this example treats a higher event rate as better. If a higher rate represented harm, the comparison would need to reflect the clinically relevant direction rather than mechanically selecting the largest value.

Monte Carlo Precision, Diagnostics, and Practical Choices

The number of posterior draws controls the simulation precision of the estimated winner fractions. For example, if an arm wins 200 of $M=1{,}000$ replicates, then $\widehat{q}_k=200/1{,}000=0.20$. Its estimated Monte Carlo standard error is

\sqrt{\frac{0.20(1-0.20)}{1{,}000}} =\sqrt{0.00016} \approx 0.01265.

Thus simulation noise in this estimated probability is about 0.01265 on the probability scale. If the same estimated probability were based on $M=10{,}000$ draws, the estimated Monte Carlo standard error would be

\sqrt{\frac{0.20(0.80)}{10{,}000}} =\sqrt{0.000016} =0.004.

A practical choice of $M$ depends on how accurately the best-arm probabilities must be estimated. Probabilities near 0.5 have larger Monte Carlo standard errors than probabilities near zero or one at the same $M$. When decisions depend on distinguishing two similar probabilities, increase $M$ and check that the estimates are stable across independent batches of draws.

Several diagnostics help detect implementation problems. Each replicate should contain exactly one winner when ties have probability zero, and the winner counts divided by $M$ should sum to one. Check that each column was sampled from the intended arm’s posterior, that the posterior parameters use the correct response coding and exposure, and that all arms use the same clinical direction. Repeating the calculation with a larger $M$ should generally reduce Monte Carlo variability, though individual estimates need not move monotonically toward their limits.

For very large $M$ or $K$, storing the entire $M$-by-$K$ matrix may use unnecessary memory. Process draws in batches: generate a block of rows, find the winners within that block, add their counts to running totals, and then discard the block. Since the total counts are sufficient for the empirical winner fractions, this gives the same calculation as storing every draw, provided the posterior draws are generated under the same model.

What This Calculation Does—and Does Not—Replace

A direct loop over arms and replicates can produce the same estimates, but it is easier to misalign arm labels or omit comparisons as $K$ grows. A vectorized row-wise maximum makes the “compare all arms in each posterior replicate” operation transparent. Vectorization is a computational organization, not a different statistical method.

Ranking posterior means answers a different question. The arm with the largest posterior mean need not have the largest probability of being best: an arm with a lower mean but a more concentrated posterior can be best less often than an arm with greater uncertainty. Likewise, selecting the arm with the largest observed response rate ignores posterior uncertainty. Thompson sampling uses posterior draws, so the relevant comparison is the sampled maximum, not either plug-in ranking.

The procedure has limits. The estimates are only as appropriate as the posterior model, response-rate definition, and independence assumptions. If posterior dependence is present, independent column draws do not represent the joint posterior and can give incorrect best-arm probabilities. The calculation also addresses “largest response rate” specifically; it does not automatically determine which treatment is best under a utility function that includes safety, cost, or other outcomes. Finally, Monte Carlo precision does not correct model misspecification or uncertainty about the clinical meaning of the endpoint.

Common mistake: “Arm B has the largest posterior mean, so its estimated best-arm probability must be one.” This is incorrect. A posterior mean is one summary of an arm’s distribution; posterior best-arm probability is the fraction of joint posterior draws in which that arm’s value exceeds every other arm’s value. Estimate it by counting row winners, not by ranking posterior means.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A matrix contains $M$ rows and $K$ columns of posterior response-rate draws. What does one row represent, and what does one column represent?
Across 500 posterior replicates, arm A is the row maximum 145 times. Calculate its estimated posterior best-arm probability.
In the preceding calculation, estimate the Monte Carlo standard error for arm A’s winner fraction.
Why is selecting the arm with the largest posterior mean not generally equivalent to estimating which arm is most likely to be best?
What basic check should hold for the winner counts across all arms when posterior draws are continuous?

Solutions. (1) A row is one joint posterior replicate containing one draw for each arm; a column is the sequence of draws for one arm. (2) The estimate is $145/500=0.290$. (3) Substitute $\widehat{q}_A=0.290$ and $M=500$: $\sqrt{0.290(0.710)/500}=\sqrt{0.0004118}\approx0.02029$. (4) The posterior mean summarizes the center of one arm’s posterior, whereas the best-arm probability counts how often that arm’s draw exceeds all other arms’ draws; posterior spread and the other arms’ distributions also matter. (5) Each replicate should have one winner, so the winner counts sum to $M$, and the estimated probabilities sum to one.