Quadrature for Nonconjugate Posterior Decision Probabilities

Why Integrate a Posterior Decision Probability?

A Bayesian decision may depend on a probability such as the posterior chance that a treatment effect is positive, exceeds a clinically relevant threshold, or is larger than another arm's effect. When the likelihood and prior do not form a conjugate pair, the posterior may not have a convenient closed form. Numerical integration provides a direct way to calculate such probabilities without replacing the model by a convenient approximation.

This tutorial uses a binary favorable outcome and a treatment effect on the log-odds scale. The calculation is illustrative: all data below are hypothetical. The method applies more broadly whenever the posterior density can be evaluated up to a normalizing constant and the decision event can be specified mathematically.

Definition: Let $\theta$ denote a model parameter, $y$ the observed data, and $A$ a region of parameter space that defines a decision event. The posterior decision probability is $$ P(\theta\in A\mid y)= \frac{\int_A L(y\mid\theta)\,\pi(\theta)\,d\theta} {\int_{\Theta}L(y\mid\theta)\,\pi(\theta)\,d\theta}, $$ where $L(y\mid\theta)$ is the likelihood, $\pi(\theta)$ is the prior density, and $\Theta$ is the parameter space. The denominator is the normalizing integral; the numerator restricts the same integrand to the decision region.

The ratio matters. A numerical approximation to the numerator alone is not a posterior probability, and the likelihood need not include constants that do not depend on $\theta$: such constants multiply both integrals and cancel. The model must nevertheless be specified correctly in both integrals.

A Nonconjugate Logistic Model

Suppose $y$ favorable outcomes are observed among $n$ participants assigned to an experimental arm. Let $\theta$ be its log odds of a favorable outcome relative to a reference log odds fixed at zero. Under the model,

Y\mid\theta\sim\operatorname{Binomial}\{n,\operatorname{logit}^{-1}(\theta)\}, \qquad \operatorname{logit}^{-1}(\theta)=\frac{1}{1+e^{-\theta}}, \qquad \theta\sim N(0,1).

The binomial likelihood is proportional, as a function of $\theta$, to $L(\theta)=p(\theta)^y\{1-p(\theta)\}^{n-y}$, where $p(\theta)=\operatorname{logit}^{-1}(\theta)$. A normal prior on a log-odds parameter and a binomial likelihood do not produce a normal posterior; this is a nonconjugate model. The assumptions for this illustration are that outcomes are conditionally independent given $\theta$, each participant has the same response probability within the arm, and the stated prior is used. The favorable-outcome coding fixes the interpretation: $\theta>0$ means a higher response probability than the reference probability of one half.

For the decision event $A=\{\theta>0\}$, define $q(\theta)=L(\theta)\phi(\theta)$, where $\phi$ is the standard normal density. Then the desired probability is the integral of $q$ over the positive half-line divided by its integral over the real line. Adaptive quadrature, such as adaptive Gauss–Kronrod integration, divides the domain into intervals and refines them where the integrand or estimated integration error requires it. The event boundary at zero is a natural place to split the integration.

Justification: Bayes' rule gives the posterior density $q(\theta)/\int_{\mathbb R}q(u)\,du$. Integrating that density over $\theta>0$ gives $$ \int_0^\infty\frac{q(\theta)}{\int_{\mathbb R}q(u)\,du}\,d\theta = \frac{\int_0^\infty q(\theta)\,d\theta}{\int_{\mathbb R}q(u)\,du}. $$ Thus the ratio of two unnormalized integrals is exactly the posterior probability; no separate posterior density formula is required.

Worked Example: Posterior Probability of a Positive Effect

Suppose, hypothetically, that 7 of 10 participants in an experimental arm have a favorable outcome. Under the model above, calculate $P(\theta>0\mid y=7)$. The binomial coefficient is omitted from the likelihood because it is constant in $\theta$ and cancels from the ratio.

For these data, $L(\theta)=p(\theta)^7\{1-p(\theta)\}^3$. The denominator and numerator are

D=\int_{-\infty}^{\infty}\phi(\theta) p(\theta)^7\{1-p(\theta)\}^3\,d\theta, \qquad N=\int_0^{\infty}\phi(\theta) p(\theta)^7\{1-p(\theta)\}^3\,d\theta.

Adaptive quadrature gives, to the displayed precision, $D\approx9.7\times10^{-4}$ and $N\approx8.4\times10^{-4}$. Therefore,

P(\theta>0\mid y=7)\approx \frac{8.4\times10^{-4}}{9.7\times10^{-4}} =0.866.

The posterior probability that the experimental arm's log odds exceed the reference log odds is about 0.866. This is a probability conditional on the model, prior, and observed data; it is not the probability that the model is true. The numerical values are rounded, so the final ratio is also approximate.

Numerical Accuracy and Integration Choices

An integration routine returns an approximation, not an exact posterior probability. A credible calculation therefore checks whether the reported value is stable under reasonable numerical refinements. Useful checks include tightening the integration tolerance, splitting the integral at the event boundary, and increasing the integration range if finite limits are used to approximate infinite tails. Results that change materially under these checks should not be treated as reliable.

For a smooth integrand over a one-dimensional domain, Gauss–Hermite quadrature is another option. With nodes $x_i$ and weights $w_i$ for the convention $\int e^{-x^2}f(x)\,dx$, a standard-normal expectation can be approximated by

\int_{-\infty}^{\infty}\phi(\theta)f(\theta)\,d\theta \approx \frac{1}{\sqrt{\pi}}\sum_i w_i f(\sqrt{2}x_i).

However, an indicator of a decision region can create a discontinuity. For $A=\{\theta>0\}$, it is often better to integrate separately over the two sides of zero, or to use adaptive quadrature, than to apply a single whole-line rule without accounting for the boundary. The denominator and numerator should use compatible accuracy settings, and a small denominator can make a ratio especially sensitive to integration error.

The same principle extends to multiple treatment effects. For a vector $\boldsymbol{\theta}$, integrate the joint likelihood times the joint prior over the decision region in the numerator and over the full parameter space in the denominator. Tensor-product quadrature applies one-dimensional rules across coordinates, but the number of grid points grows rapidly with dimension. Adaptive or sparse-grid methods may help in some settings, yet no method removes the need to check numerical stability and the behavior of the integrand near decision boundaries.

Worked Example: Comparing Two Identically Modeled Arms

Consider two hypothetical experimental arms with independent effects $\theta_1,\theta_2$, each having the $N(0,1)$ prior. Suppose each arm has 6 favorable outcomes among 10 participants, and outcomes in each arm follow the logistic binomial model above. The decision event is that arm 1 has the larger log-odds effect: $A=\{\theta_1>\theta_2\}$.

Let $q(\theta)=\phi(\theta)p(\theta)^6\{1-p(\theta)\}^4$, and let $D=\int_{\mathbb R}q(\theta)\,d\theta$. The posterior factors into two identical one-dimensional densities, so its normalizing integral is $D^2$. The numerator is

N=\int_{\theta_1>\theta_2}q(\theta_1)q(\theta_2)\,d\theta_1d\theta_2.

Swapping the two coordinates maps the region $\theta_1>\theta_2$ onto $\theta_2>\theta_1$ and leaves the integrand unchanged. The diagonal $\theta_1=\theta_2$ has probability zero under this continuous posterior. The two regions therefore have equal integrals, and together they contain the full integral $D^2$. Hence $N=D^2/2$, giving

P(\theta_1>\theta_2\mid y_1=y_2=6)=\frac{N}{D^2}=\frac{1}{2}=0.500.

A two-dimensional quadrature implementation can integrate over one triangular region and confirm the result, or integrate over both regions with the diagonal handled consistently. This example is a useful computational check: identical data, priors, and likelihoods should not produce a systematic preference for either arm. It does not imply that the effects are equal; it says each is equally likely to be larger under this symmetric posterior.

Quadrature Compared With Simulation

Quadrature evaluates an approximation to an integral directly. Simulation instead approximates a posterior probability by the fraction of posterior draws that fall in the decision region. The earlier material on Monte Carlo error in posterior probabilities gives the sampling-error framework for interpreting that fraction. Quadrature has no Monte Carlo error, but it does have numerical integration error; simulation has Monte Carlo error and may also have sampling or convergence problems, particularly for correlated draws.

For a low-dimensional, smooth posterior and a clearly defined event, quadrature can be fast, deterministic, and highly accurate. Simulation is often more flexible when the parameter dimension is large or the posterior region is irregular. A simulation estimate can also provide an independent check on a quadrature result. Neither method is automatically trustworthy: quadrature requires convergence and tail checks, while simulation requires adequate exploration of the posterior and an appropriate estimate of effective sample size when draws are dependent.

Worked Example: Checking Quadrature Against Posterior Simulation

Return to the hypothetical 7-of-10 example. Suppose a simulation procedure produces 100,000 effectively independent draws from the specified posterior, of which 86,600 have $\theta>0$. The simulation estimate is

\widehat{P}(\theta>0\mid y=7)=\frac{86{,}600}{100{,}000}=0.866.

Treating these as independent Bernoulli indicators, the estimated Monte Carlo standard error is

SE_{\mathrm{MC}} = \sqrt{\frac{0.866(1-0.866)}{100{,}000}} = \sqrt{0.000001161} \approx0.00108.

The quadrature value, approximately 0.866, agrees with the simulation estimate to the reported precision. The simulation's approximate 95% Monte Carlo range is $0.866\pm1.96(0.00108)$, or about 0.864 to 0.868. This range describes simulation uncertainty in estimating the posterior probability; it is not a credible interval for $\theta$. If the retained draws are autocorrelated, the raw draw count is not the effective independent sample size, and the Monte Carlo error must account for that dependence.

Diagnostics, Limitations, and Common Pitfalls

A probability calculation can look plausible while still being numerically wrong. Check that the denominator is positive and finite, that the numerator is between zero and the denominator, and that the resulting ratio lies in $[0,1]$. Refine the quadrature rule or tolerance and compare the result; inspect tail contributions or expand any finite integration limits; and split the domain at discontinuities such as an event boundary. In multiple dimensions, compare increasingly accurate grids or alternative integration strategies when feasible.

Using a numerator from a different integrand. The numerator must use the same likelihood and prior as the denominator, restricted only to the decision region.
Ignoring a decision boundary. A numerical rule can poorly approximate a discontinuous indicator if the boundary crosses its grid. Split the integration region or use a method designed to handle that boundary.
Confusing numerical accuracy with model adequacy. A stable integral is accurate for the specified model and prior; numerical precision does not establish that those assumptions are appropriate.
Expecting ordinary tensor quadrature to scale indefinitely. Grid size grows quickly as the number of parameters increases. Simulation or other integration strategies may be preferable in higher dimensions.

A common mistake is to report a numerator of $6.2\times10^{-4}$ with a denominator of $9.7\times10^{-4}$ for the 7-of-10 example and conclude that the posterior probability is about 0.64. The ratio arithmetic is correct for those two numbers, but the numerator does not match the specified integral. Re-evaluating that integral gives approximately $8.4\times10^{-4}$, so the posterior probability is about 0.866. Always validate the integral values themselves, not only their division.

Quadrature is most useful when a decision probability is low-dimensional enough to integrate accurately and reproducibly. Its central benefit is a direct, deterministic calculation; its central limitation is computational growth with dimension and sensitivity to poorly handled tails or boundaries. Simulation is a flexible alternative, and comparing the two can reveal errors that either approach alone may not make apparent.

Check Your Understanding

Work through each question, then compare with the solutions below.

In the posterior probability formula, what distinguishes the numerator from the denominator?
If a numerical calculation gives a denominator of 0.012 and a numerator of 0.009, what posterior probability does it imply?
A simulation estimate is 0.42 from 2,500 effectively independent draws. Calculate its estimated Monte Carlo standard error.
Why should a two-arm quadrature calculation with identical priors, data, and likelihoods give a probability of 0.5 that arm 1 exceeds arm 2?

Solutions. (1) The numerator integrates the likelihood times prior only over the decision region; the denominator integrates the same product over the full parameter space. (2) The ratio is $0.009/0.012=0.75$. (3) The standard error is $\sqrt{0.42(0.58)/2500}=\sqrt{0.00009744}\approx0.00987$. (4) Swapping the two arm parameters leaves the joint posterior unchanged and swaps the two strict-ordering regions. The diagonal has probability zero, so each region has probability one half.