Ratio of Uniforms: The Theorem

Theorem (Kinderman-Monahan, 1977). Let \(f(x) \geq 0\) be a density (possibly unnormalized). Define region: \[\mathcal{R}_f = \left\{(u,v): 0 \leq u \leq \sqrt{f(v/u)}\right\}\]

If \((U,V) \sim \Uniform(\mathcal{R}_f)\), then \(X = V/U\) has density proportional to \(f\).

In words: scatter points uniformly on a region in the \((u,v)\) plane defined by \(f\), then read off \(X = V/U\). The transformation from uniform samples to the target happens through the geometry of the region.

Ratio of Uniforms: Derivation

Change of variables. \((u,v) \to (u, x)\) where \(x = v/u\), so \(v = ux\).

Jacobian. \(\dfrac{\partial(u,v)}{\partial(u,x)} = \begin{pmatrix} 1 & 0 \\ x & u \end{pmatrix}\), so \(|J| = u\).

Joint density of \((U,X)\). Since \((U,V)\) is uniform on \(\mathcal{R}_f\): \[p(u,x) = \frac{u}{|\mathcal{R}_f|} \cdot \mathbf{1}\left[0 \leq u \leq \sqrt{f(x)}\right]\]

Marginal of \(X\): \[p(x) = \int_0^{\sqrt{f(x)}} \frac{u}{|\mathcal{R}_f|} \, du = \frac{1}{|\mathcal{R}_f|} \cdot \frac{f(x)}{2} \propto f(x) \qquad \square\]

Finding the Bounding Rectangle

The Challenge. We need to sample uniformly from \(\mathcal{R}_f\). Solution: find a rectangle containing \(\mathcal{R}_f\).

Ratio of Uniforms: Examples

ARS: The Problem

Recall: Accept-Reject Challenges. For complex distributions, finding a good proposal \(g\) is hard.

• Too loose: poor efficiency
• Too tight: difficult to sample from
• Fixed envelope: can’t improve

Log-Concavity

Definition. A density \(f\) is log-concave if \(\log f\) is concave: \[\log f(\lambda x + (1-\lambda)y) \geq \lambda \log f(x) + (1-\lambda)\log f(y)\]

ARS: Initial Envelope

• 4 initial points: Tangents create loose upper hull
• Shaded region: Gap between upper and lower bounds
• Efficiency: Proportional to gap area (smaller is better)

ARS: After Adaptation

• 7 support points: Added 3 rejected samples
• Tighter envelope: Shaded region much smaller
• Higher efficiency: proposals per accepted draw fall from 1.08 to 1.04; each rejection improved the envelope for good

ARS: Advantages and Extensions

Monte Carlo Integration: A Recap

The Problem. Evaluate: \(I = \int_{\mathcal{X}} h(x) f(x) dx = \E_f[h(X)]\)

Basic Monte Carlo. Sample \(X_1, \ldots, X_n \iid f\), estimate \(\hat{I}_n = \frac{1}{n}\sum_{i=1}^n h(X_i)\).

• \(\E[\hat{I}_n] = I\) (unbiased)
• \(\Var[\hat{I}_n] = \frac{\Var_f[h(X)]}{n}\)
• \(\hat{I}_n \xrightarrow{a.s.} I\) (SLLN)
• \(\sqrt{n}(\hat{I}_n - I) \xrightarrow{d} \Normal(0, \Var_f[h(X)])\)

Recall the \(n^{-1/2}\) rate from week 1? It is dimension-free: the same rate in one dimension and in a million. Everything that follows is about pushing the constant in front of \(n^{-1/2}\) down.

Error Rate. Standard error: \(O(n^{-1/2})\), independent of dimension. Compare: Riemann sums \(O(n^{-1/d})\).

Importance Sampling: Unbiasedness and Variance

Unbiasedness. \[\E_g[\hat{I}_{IS}] = \E_g\left[h(X)\frac{f(X)}{g(X)}\right] = \int h(x)\frac{f(x)}{g(x)}g(x)dx = \int h(x)f(x)dx = I\]

Variance. \[\Var_g[\hat{I}_{IS}] = \frac{1}{n}\Var_g\left[h(X)\frac{f(X)}{g(X)}\right] = \frac{1}{n}\left(\int \frac{h(x)^2 f(x)^2}{g(x)} dx - I^2\right)\]

The variance depends on the match between \(g\) and \(|h|f\): where \(g\) underweights regions that matter for \(hf\), the importance weights blow up and the variance grows.

Finite-variance condition. The estimator has finite variance when \(\int \frac{h(x)^2 f(x)^2}{g(x)} dx < \infty\).

Importance Sampling: Optimal Proposal

Optimal Proposal. Minimum variance when: \[g^*(x) \propto |h(x)|f(x)\]

Zero variance if \(h \geq 0\). But \(g^*\) requires knowing \(\int |h(x)|f(x)dx\), which is the integral you are trying to compute.

When Importance Sampling Fails

Left: estimating \(P(X > 3)\) for a standard normal. A \(\Normal(3,1)\) proposal concentrates draws where the integrand lives and beats plain MC by orders of magnitude. Right: estimating \(\E[X^2]\) for \(X \sim t_3\). A \(\Normal(0,1)\) proposal has lighter tails, so the weights are unbounded and the running estimate never settles. A \(t_1\) proposal has heavier tails and converges.

Control Variates

The Idea. Use a correlated variable \(Y\) with known \(\E[Y] = \mu_Y\): \[\hat{I}_{CV} = \hat{I} - c(\bar{Y} - \mu_Y)\]

Optimal \(c\). \[c^* = \frac{\Cov(h(X), Y)}{\Var(Y)}\]

Variance reduction: \[\frac{\Var[\hat{I}_{CV}]}{\Var[\hat{I}]} = 1 - \rho^2_{h(X),Y}\]

The optimal coefficient \(c^*\) is the slope of the regression of \(h(X)\) on \(Y\). The variance reduction equals the \(R^2\) of that regression: the more \(Y\) explains \(h(X)\), the more variance you remove.

Example. Estimating \(\E[e^X]\) where \(X \sim \Normal(0,1)\): use \(Y = X\) as control, \(\E[Y] = 0\) known. High correlation.

Control Variate Adjustment

Antithetic Variables

The Idea. Generate pairs \((X, X')\) with:

• Same marginal distribution
• Negative correlation

Implementation. For \(U \sim \Uniform(0,1)\):

• Use both \(U\) and \(1-U\)
• Transform both via \(F^{-1}\)

Estimate: \(\hat{I}_{AV} = \frac{1}{2n}\sum_{i=1}^n [h(X_i) + h(X'_i)]\)

Antithetic Variables: Detailed Analysis

Variance Reduction Formula. For antithetic pairs \((X, X')\) with \(X' = F^{-1}(1-U)\) when \(X = F^{-1}(U)\): \[\Var[\hat{I}_{AV}] = \frac{1}{2n}[\Var(h(X)) + \Cov(h(X), h(X'))]\]

Reduction factor: \(\dfrac{\Var[\hat{I}_{AV}]}{\Var[\hat{I}_{MC}]} = \dfrac{1}{2}[1 + \rho_{h(X),h(X')}]\)

The factor compares \(n\) antithetic pairs with \(n\) independent draws. A pair costs two evaluations of \(h\), so it has to beat \(1/2\) to pay for itself: it does when \(h\) is monotone (\(\rho < 0\)), and it loses when \(\rho > 0\), as for a U-shaped or oscillating \(h\).

Antithetic Variables: Example

Example: \(\E[e^X]\), \(X \sim \Uniform(0,1)\).

• \(X' = 1 - X\) (perfect negative correlation)
• \(\Cov(e^X, e^{1-X}) = e - (e-1)^2 \approx -0.234\)
• \(\Var(e^X) = (e^2-1)/2 - (e-1)^2 \approx 0.242\)
• \(\rho_{h(X),h(X')} \approx -0.968\)
• Reduction factor: \(\frac{1}{2}(1 - 0.968) = 0.016\)

\(e^X\) is nearly linear on \([0, 1]\), so this is close to the best case; the next slide gives the typical range.

Antithetic Variables: Practical Gains

When Antithetic Works Best.

• \(h\) is monotonic: \(\rho < 0\)
• Linear \(h\) with symmetric \(f\): perfect
  • \(\rho = -1\)
  • Reduction factor = 0
  • Zero variance.

Typical reduction factors:

• Monotonic \(h\): 0.1–0.3
• Non-monotonic \(h\): 0.4–0.6
• Worst case: \(> 0.5\) (harmful)

Stratified Sampling: Detailed Analysis

Variance Reduction Formula. With \(k\) strata, proportions \(p_j = P(X \in S_j)\), and \(n_j\) samples from stratum \(j\): \[\Var[\hat{I}_{SS}] = \sum_{j=1}^k p_j^2 \frac{\sigma_j^2}{n_j}\]

Four Variance-Reduction Ideas, Four Roles

Variance Reduction Comparison

All four methods spend the same 1000 evaluations of \(h\), so the antithetic ratio shown is \(1 + \rho = 0.032\), twice the per-pair factor \(\tfrac12(1+\rho) = 0.016\) quoted earlier; the control variate \(Y = U\) (\(\rho = 0.992\)) and ten proportional strata do about as well on this near-linear integrand.

The Ugly Truth

What Textbooks Don’t Tell You. Those beautiful variance reduction formulas assume we know things we don’t.

The Setup Costs.

• Stratified: Need pilot study to estimate \(p_j\), \(\sigma_j\)
• Control Variates: Must know \(\E[Y]\) exactly
• Importance Sampling: Finding good proposal is hard
• Antithetic: Only works if \(h\) is monotonic

Reality Check. 10 strata, 1000 pilot samples, 10000 main samples:

• 10% overhead just for setup
• Need more than 10% variance reduction to break even
• Often don’t get it.

When Classical Methods Fail

Stratified Sampling Dies In:

• High dimensions (curse of dimensionality)
• Complex distributions (can’t compute \(P(X \in S_j)\))
• Black-box simulators (no access to internals)

Importance Sampling Dies In:

• High dimensions (weight degeneracy)
• Tail mismatch: proposal lighter than the target (infinite variance)
• Bad proposal choice (most of them)

The Brutal Truth. In practice, just using more samples often beats fancy variance reduction.

PSIS

Pareto Smoothed Importance Sampling. The one method that actually works in practice.

The Key Idea. Extreme importance weights follow a Pareto distribution.

• Fit generalized Pareto to largest 20% of weights
• Replace extreme weights with smoothed values
• Get automatic reliability diagnostic (\(\hat{k}\))

The \(\hat{k}\) Diagnostic.

• \(\hat{k} < 0.5\): Excellent, trust results
• \(0.5 \leq \hat{k} < 0.7\): OK, use with caution
• \(\hat{k} \geq 0.7\): Bad, don’t trust, recompute

Why PSIS Dominates

Where You’ll See It.

• Stan (loo package)
• ArviZ (Python/R)
• Any modern Bayesian software

The Killer App: LOO-CV from One MCMC Run. Without PSIS: Run MCMC \(n\) times (impossible). With PSIS: Run MCMC once, get LOO-CV for free.

But more on this when we formally enter the realms of Bayesian computing.

What to Carry Out of the Room

• Ratio of uniforms replaces Accept-Reject’s proposal \(g\) and bound \(M\) with bounds on a 2D region, and handles unnormalized densities automatically.
• ARS builds and tightens its own envelope as it runs, for any log-concave target.
• Importance sampling is the centerpiece: keep every draw, weight by \(f/g\), and the variance of the weighted average depends entirely on how well \(g\) matches \(|h|f\).
• The four variance-reduction techniques each exploit a different structural feature and each can backfire when the feature is absent.

The week’s punchline: once you move from generating samples to estimating integrals, the design problem shifts from “which draws do I keep?” to “which draws do I weight, and how?” Importance sampling is the bridge between the sampling methods of weeks 2–4 and the Bayesian computing machinery from week 10 on, and the \(\hat{k}\) diagnostic tells you whether the bridge is holding.

Next week: the bootstrap. A different answer to “how do I estimate variability?”, one that replaces distributional assumptions with the empirical distribution.