Checking Posterior Best Probabilities by Numerical Integration

Why Check a Posterior Best Probability?

In Thompson sampling, a treatment’s assignment probability can depend on its posterior probability of being best. For two treatments, this probability is conceptually simple, but its calculation may require integrating over a posterior distribution. A numerical result can look precise while still being wrong because of a coding error, an unsuitable integration range, or a missed feature of the posterior density. Checking the integral against a separate Monte Carlo estimate is a practical way to detect such problems.

The goal here is to calculate the probability for two treatments using numerical quadrature: approximating an integral from values of its integrand at selected points. Monte Carlo uses a different approach: repeatedly drawing parameter values from the posterior and counting how often one treatment is best. The methods have different sources of numerical error, which makes their comparison useful.

From Posterior Probability to an Integral

Definition: Let $\theta_1$ and $\theta_2$ be the treatment parameters on a utility scale where larger values are better, and let $D$ denote the observed trial history. The posterior best probability for treatment 1 is $q_1=P(\theta_1>\theta_2\mid D)$. If the two posterior distributions are independent, with densities $f_1$ and $f_2$ and cumulative distribution function $F_2$ for $\theta_2$, then $$ q_1=\int_{-\infty}^{\infty} f_1(x)F_2(x)\,dx. $$ For parameters restricted to a range, the integral is taken over that range.

The independence condition is important. It is satisfied, for example, when the model and prior give independent treatment parameters and the observations update each treatment separately. If the posterior parameters are dependent, the product of a marginal density and a marginal cumulative distribution function is not generally the correct integrand. In that case, the corresponding expression is $q_1=\int f_1(x)P(\theta_2<x\mid\theta_1=x,D)\,dx$.

For independent, continuous posteriors, the probability of a tie is zero, so $q_2=P(\theta_2>\theta_1\mid D)=1-q_1$. If the posterior can place positive probability on a tie, the strict inequality and tie rule must be handled explicitly; treating ties as impossible would change the answer. We use the larger-is-better convention throughout. If smaller outcomes are preferable, define a utility that reverses the direction before calculating the best probability.

Why the Integral Has This Form

Condition on the value of $\theta_1$. If $\theta_1=x$, treatment 1 is best when $\theta_2<x$, which has posterior probability $F_2(x)$. Averaging this conditional probability over the posterior density of $\theta_1$ gives the result.

\begin{aligned} P(\theta_1>\theta_2\mid D) &=\int P(\theta_1>\theta_2\mid\theta_1=x,D) f_1(x)\,dx\\ &=\int P(\theta_2<x\mid D) f_1(x)\,dx\\ &=\int F_2(x)f_1(x)\,dx. \end{aligned}

The second equality uses posterior independence. The derivation also explains why numerical quadrature is a natural tool: it evaluates a one-dimensional integral even when no convenient closed-form probability is available.

Numerical Quadrature in Practice

A quadrature rule approximates an integral using a weighted sum of integrand values. For a smooth function $g$ on $[a,b]$, composite Simpson quadrature divides the interval into an even number $n$ of equal subintervals, with width $h=(b-a)/n$. Writing $x_j=a+jh$, its approximation is

\int_a^b g(x)\,dx\ \approx\ \frac{h}{3}\left[g(x_0)+g(x_n) +4\sum_{\substack{j=1\\ j\ \mathrm{odd}}}^{n-1}g(x_j) +2\sum_{\substack{j=2\\ j\ \mathrm{even}}}^{n-2}g(x_j)\right].

The rule is exact for polynomials of degree at most three, but that fact does not guarantee accuracy for every posterior integrand. For smooth densities on unbounded ranges, adaptive Gauss–Kronrod quadrature is often convenient; for finite intervals and a transparent hand calculation, composite Simpson quadrature is useful. Always check whether the integration range covers the posterior mass and whether the integrand is smooth enough for the chosen rule.

Worked Example: Beta Posteriors and Simpson Quadrature

Suppose two hypothetical treatments have independent Bernoulli outcomes, with success probability $\theta_i$ for treatment $i$. Use independent $\mathrm{Beta}(1,1)$ priors. The hypothetical trial history contains one success and no failures on treatment 1, and no successes and one failure on treatment 2. The beta-binomial update gives $\theta_1\mid D\sim\mathrm{Beta}(2,1)$ and $\theta_2\mid D\sim\mathrm{Beta}(1,2)$. The Bernoulli model and independent priors give independent posteriors, and larger success probability is better.

On $[0,1]$, the posterior density of $\theta_1$ is $f_1(x)=2x$, and the posterior cumulative distribution function of $\theta_2$ is $F_2(x)=1-(1-x)^2=2x-x^2$. Thus $g(x)=f_1(x)F_2(x)=4x^2-2x^3$, and $q_1=\int_0^1g(x)\,dx$.

Use Simpson quadrature with $n=4$, so $h=(1-0)/4=0.25$. The nodes are $0,0.25,0.50,0.75,1$; substituting them gives $g(x_j)=0,0.21875,0.75000,1.40625,2.00000$. Therefore,

q_1\approx \frac{0.25}{3}\left[0+2+4(0.21875+1.40625)+2(0.75000)\right] =\frac{10}{12}=0.8333.

In fact, integrating the cubic exactly gives $q_1=5/6=0.8333$; Simpson’s rule is exact here because the integrand has degree three. For an independent computational check, suppose a hypothetical Monte Carlo run draws 10,000 posterior pairs and finds treatment 1 larger in 8,331 pairs. The estimate is $8{,}331/10{,}000=0.8331$, with Monte Carlo standard error approximately $\sqrt{0.8331(1-0.8331)/10{,}000}=0.00373$. The difference from the quadrature value, $0.00023$, is small relative to this sampling error.

Worked Example: Normal Posteriors on an Unbounded Range

Suppose hypothetical independent normal posteriors are $\theta_1\mid D\sim N(0.60,0.04)$ and $\theta_2\mid D\sim N(0.50,0.09)$. The second argument is the variance, so the posterior standard deviations are 0.20 and 0.30. Assume the model is appropriate and the two posteriors are independent. The best probability is the integral of $f_1(x)F_2(x)$ over the real line.

For numerical integration, set $x=0.60+0.20z$. The transformed integral is $\int_{-\infty}^{\infty}\phi(z)\Phi((0.10+0.20z)/0.30)\,dz$, where $\phi$ and $\Phi$ are the standard normal density and cumulative distribution function. An adaptive quadrature calculation over $[-6,6]$ gives approximately $0.6092$; the omitted standard-normal tails have total probability less than $2.0\times10^{-9}$.

A useful check is available because the difference of independent normal variables is normal: $\theta_1-\theta_2$ has mean $0.60-0.50=0.10$ and variance $0.04+0.09=0.13$. Consequently,

P(\theta_1>\theta_2\mid D) =P(\theta_1-\theta_2>0\mid D) =\Phi\left(\frac{0.10}{\sqrt{0.13}}\right) =\Phi(0.2774)\approx 0.6092.

The numerical integral agrees with this calculation. As a separate hypothetical Monte Carlo check, suppose 60,941 of 100,000 independent posterior draws have $\theta_1>\theta_2$. The estimate is $0.60941$, and its estimated standard error is $\sqrt{0.60941(1-0.60941)/100{,}000}\approx0.00154$. This agreement is reassuring, but it does not by itself prove that either implementation is correct.

Worked Example: A Change of Variable at a Singular Endpoint

Consider hypothetical independent posteriors $\theta_1\mid D\sim\mathrm{Beta}(1/2,1)$ and $\theta_2\mid D\sim\mathrm{Beta}(1,1)$, with larger values better. These are valid beta distributions on $[0,1]$. Here $f_1(x)=\tfrac12x^{-1/2}$, while $F_2(x)=x$. The integrand $g(x)=\tfrac12\sqrt{x}$ is finite, but its derivative becomes unbounded as $x$ approaches zero, which can make a uniform-grid rule converge less smoothly.

The direct integral is $q_1=\int_0^1\tfrac12\sqrt{x}\,dx=1/3$. To improve the numerical behavior, set $x=t^2$, so $dx=2t\,dt$. The transformed integral becomes $\int_0^1 t^2\,dt$. With Simpson quadrature using $n=2$ and $h=0.5$, the node values are $0,0.25,1$, giving

q_1\approx \frac{0.5}{3}\left[0+1+4(0.25)\right] =\frac{1}{3}=0.3333.

The substitution has made the integrand smooth and, in this case, the rule is exact. For a hypothetical Monte Carlo check, suppose 3,336 of 10,000 posterior pairs favor treatment 1. The estimate is $0.3336$, with standard error $\sqrt{0.3336(1-0.3336)/10{,}000}\approx0.00471$, consistent with the quadrature result.

Diagnostics, Alternatives, and Limitations

A good numerical check looks for agreement that is plausible given the error of each method, not identical digits. Monte Carlo estimates fluctuate between runs; quadrature estimates can have deterministic error from the grid, integration range, or implementation. Useful checks include:

Refine the quadrature. Increase the number of subintervals or tighten the adaptive integration tolerance. Stable results across refinements are more credible than a result from one setting.
Check the range and tails. For an unbounded posterior, use a justified transformation or verify that omitted tail probability is negligible at the reported precision.
Check the posterior inputs. Densities should integrate to one, and the CDF should be nondecreasing with values between zero and one. Verify that the model and trial history produce the stated posterior parameters.
Check the probability constraints. The result must lie between zero and one. With independent continuous posteriors, separately calculated probabilities for each arm should sum to approximately one; a substantial discrepancy signals error or inadequate integration accuracy.
Compare at the right precision. For a Monte Carlo estimate $\hat q$ from $M$ independent posterior pairs, its estimated standard error is $\sqrt{\hat q(1-\hat q)/M}$. A difference of a few such standard errors deserves investigation, while small differences are expected.

Quadrature is deterministic for fixed settings and can be very accurate for a smooth one-dimensional integrand. It can be preferable to Monte Carlo when high precision is needed and the integral is manageable. Monte Carlo is straightforward to implement from posterior draws and extends naturally to settings where integration is difficult, but its error typically decreases only at rate $1/\sqrt{M}$. For two independent arms with a known closed-form comparison, that analytic calculation may be the simplest benchmark. These approaches are complementary: use an available analytic result to validate quadrature, or use Monte Carlo as a computational cross-check when no closed form is available.

Numerical integration is not a substitute for a correct posterior model. It cannot repair an incorrect likelihood, a mistaken treatment coding, an inappropriate utility scale, or an independence assumption that does not hold. Nor does a posterior best probability alone establish that the treatment has clinically meaningful benefit; it answers which parameter is larger under the specified posterior and utility definition.

Common pitfall: Treating agreement between quadrature and Monte Carlo as proof that the posterior calculation is correct. Both methods can use the same incorrectly specified posterior or compare the wrong direction of benefit. Check the model, treatment coding, parameter orientation, and posterior dependence separately; then use agreement as evidence about the numerical calculation.

Check Your Understanding

Answer each question, then compare with the brief solutions.

For independent continuous posteriors, explain why $P(\theta_2>\theta_1\mid D)=1-P(\theta_1>\theta_2\mid D)$.
A quadrature estimate is 0.72, and a Monte Carlo estimate is 0.70 from 400 independent posterior pairs. Is this disagreement necessarily concerning? Use the Monte Carlo standard error to explain.
Why is $f_1(x)F_2(x)$ generally not the correct integrand when the posterior parameters are dependent?
What diagnostic is especially important when integrating a normal-posterior expression over a finite interval?

Solutions. (1) Continuous independent posteriors have zero probability of a tie, so one of the two strict orderings must hold. (2) The Monte Carlo standard error is approximately $\sqrt{0.70(0.30)/400}=0.0229$; the difference $0.02$ is less than one standard error and is not, by itself, concerning. (3) Dependence means the conditional probability of $\theta_2<x$ can depend on $\theta_1=x$; use $f_1(x)P(\theta_2<x\mid\theta_1=x,D)$. (4) Check that omitted tail mass is negligible, for example by enlarging the interval and confirming that the result is stable.