Probability assigns events a number in \([0,1]\):
Example: dice roll, \(P(\{4\})=1/6\).
For continuous variables, we require an integral for it to make sense:
\[ P(3 \le x \le 10) = \int_3^{10}p(x)\,dx =0.95. \]
Probability = area under the probability density
\[ (\Omega,\mathcal F,P) \]
\(\Omega\): outcome / parameter space
\(\mathcal F\): measurable events \(A\subseteq\Omega\)
\(P\): probability measure
\[ P(A)\ge0,\qquad P(\Omega)=1. \]
For disjoint \(A_i\):
\[ P\!\left(\bigcup_i A_i\right)=\sum_iP(A_i). \]
\[ \Omega=\mathbb R^n. \]
We integrate with respect to Lebesgue measure on \(\mathbb R^n\), written \(dx\).
From now on, densities are with respect to this measure:
\[ p(x)=\frac{dP}{dx}(x), \qquad P(A)=\int_A p(x)\,dx. \]
Define a Bayesian model: prior + likelihood \(\rightarrow\) joint density.
Prior choice: \(p(x\mid M)\) \(\propto\) \(\sigma^{-1}\)
Likelihood for \(D\): \(p(D\mid x,M)\) \(=\) \((2\pi\sigma^2)^{-3/2}\prod_{i=1}^3\exp\!\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right]\)
Joint density:
| \(p(D,x\mid M)\) | \(=\) | \(p(D\mid x,M)\) \(p(x\mid M)\) |
| \(p(D,x\mid M)\) | \(\propto\) | \(\sigma^{-3}\prod_{i=1}^3\exp\!\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right]\) \(\sigma^{-1}\) |
Parameter space:
\[ x=(\mu,\sigma)\in\mathbb R\times\mathbb R_+. \]
Factorization:
\[ p(D,x\mid M) = p(D\mid x,M)p(x\mid M). \]
Likelihood contribution:
\[ p(D\mid x,M)\propto_x \sigma^{-3} \prod_{i=1}^3 \exp\!\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right]. \]
Prior contribution:
\[ p(x\mid M)\propto\sigma^{-1}. \]
Therefore:
\[ p(D,x\mid M) \propto_x \sigma^{-4} \prod_{i=1}^3 \exp\!\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right]. \]
Bayes, “An essay towards solving a problem in the doctrine of chances”, Philos. Trans. R. Soc. Lond., 1763.
Same joint density:
\[ p(D,x\mid M). \]
Condition on \(x\): likelihood.
Condition on \(D\): posterior.
\(p(D\mid x,M)\) \(=\) \(p(D,x\mid M)\) \(p(x\mid M)\)
But this is the term we defined before.
\(p(x\mid D,M)\) \(=\) \(p(D\mid x,M)\) \(p(x\mid M)\) \(p(D\mid M)\)
Bayes
\(p(D\mid M)\) \(=\) \(\int\) \(p(D\mid x,M)\) \(p(x\mid M)\) \(dx\)
Joint:
\[ p(D,x\mid M) = p(D\mid x,M)p(x\mid M). \]
Likelihood:
\[ p(D\mid x,M) = \frac{p(D,x\mid M)} {\int p(D,x\mid M)\,dD}. \]
Parameter conditional:
\[ p(x\mid D) = \frac{p(D,x\mid M)} {\int p(D,x\mid M)\,dx}. \]
For this model:
\[ dx=d\mu\,d\sigma. \]
Bayes, “An essay towards solving a problem in the doctrine of chances”, Philos. Trans. R. Soc. Lond., 1763.
Example: is variable \(x_3\) larger than \(0\)?
\[ P(x_3>0\mid D) = \int_{\{x:x_3>0\}} p(x\mid D)\,dx. \]
Many posterior questions become calculations of the form:
\[ \mathbb E_p[f] = \int_\Omega f(x)\,p(x\mid D)\,dx. \]
\(f(x) = x\) calculate the mean of the variable \(x\).
Variable/parameter:
\[ x\in\Omega. \]
Observable:
\[ f:\Omega\to\mathbb R. \]
Second moment:
\[ \mathbb E_p[x_j^2] = \int x_j^2\,p(x\mid D)\,dx. \]
Interval probability:
\[ f(x)=\mathbf 1_A(x), \qquad P(x\in A\mid D)=\mathbb E_p[f]. \]
Credible interval \(C_{1-\alpha}\):
\[ F_j(t) = P(x_j\le t\mid D). \]
\[ \begin{aligned} F_j(a)&=\alpha/2,\\ F_j(b)&=1-\alpha/2 \end{aligned} \quad\text{solve for }a,b. \]
\[ C_{1-\alpha}=[a,b]. \]
Easiest approach: Riemann sum.
Replace area by small rectangles:
\[
\int f(x)p(x)\,dx
\approx
\sum_i f(x_i)p(x_i)\Delta x.
\] 
More intelligent nodes and weights: Gaussian quadrature.
Excellent in low dimension when the integrand is polynomial-like.
Integrand:
\[ g(x)=f(x)p(x\mid D). \]
Riemann:
\[ \int g(x)\,dx \approx \sum_i g(x_i)\Delta x. \]
Gaussian quadrature:
\[ \int_a^b g(x)\rho(x)\,dx \approx Q_n(g) := \sum_{i=1}^n \lambda_i g(x_i). \]
Choose \(x_i\) as roots of the degree-\(n\) polynomial orthogonal w.r.t. \(\rho\).
Optimal in polynomial class: orthogonality cancels the interpolation-error moments, so \(Q_n\) matches \(2n\) moments:
\[ Q_n(g)=\int_a^b g(x)\rho(x)\,dx \quad\text{for }\operatorname{degree}(g)\le 2n-1. \]
Riemann, “Ueber die Darstellbarkeit einer Function durch eine trigonometrische Reihe”, Abh. K. Ges. Wiss. Gottingen, 1868. Gauss, “Methodus nova integralium valores per approximationem inveniendi”, Comment. Soc. Regiae Sci. Gotting. Recent., 1814.
Better deterministic idea: refine intervals/cells only when needed.
Use a local error estimate:
On a cell \(I\):
\[ Q_I(g)\approx\int_I g(x)\,dx. \]
Embedded estimate:
\[ \widehat e_I = |Q_I^{\mathrm{high}}(g)-Q_I^{\mathrm{low}}(g)|. \]
Refine if:
\[ \widehat e_I>\tau_I. \]
Examples:
Still partitions ambient space.
Kuncir, “Algorithm 103: Simpson’s rule integrator”, Commun. ACM, 1962. Gander et al., “Adaptive Quadrature - Revisited”, BIT Numer. Math., 2000.
Mesh depends on several choices/parameters.
Simpler/more elegant: random delta peaks.
\[ p(x)\,dx \approx \sum_i w_i\delta_{x_i}(dx). \]
Draw from an easier proposal:
\[ x_i\sim q. \]
Correct by density ratios:
\[ w_i = \frac{p(x_i)/q(x_i)} {\sum_j p(x_j)/q(x_j)}. \]
Failure: proposal misses support or tails.
For \(N\) proposals:
\[ x_i\sim q, \qquad r_i=\frac{p(x_i)}{q(x_i)}. \]
Normalized weights:
\[ \bar w_i = \frac{r_i}{\sum_j r_j}. \]
Useful fraction:
\[ \eta_{\mathrm{IS},N} = \frac{N_{\mathrm{eff}}}{N} \approx \frac{(\sum_i r_i)^2} {N\sum_i r_i^2}. \]
Large \(N\):
\[ \eta_{\mathrm{IS}} \to \frac{Z^2}{\mathbb E_q r^2}, \qquad Z=\mathbb E_q r. \]
If \(p\) normalized: \(\mathbb E_q r=1\).
Kahn, “Random sampling (Monte Carlo) techniques in neutron attenuation problems–I”, Nucleonics, 1950. Kish, “Survey Sampling”, Wiley, 1965.
Weighted samples: correct, but not a real sampler from \(p\).
Unweighted samples by accept/reject:
\[ x\sim q, \qquad u\sim \mathrm{Unif}(0,1), \qquad u\le \frac{p(x)}{c\,q(x)}. \]
Need valid probabilities:
\[ c\,q(x)\ge p(x) \quad\text{for all }x. \]
Retained fraction \(\approx Z/c\); cost \(\approx c/Z\) proposals per sample.
For \(N\) proposals:
\[ x_i\sim q, \qquad r_i=\frac{p(x_i)}{q(x_i)}, \qquad Z=\mathbb E_q r. \]
Envelope:
\[ r_i\le c \quad\text{for all proposals.} \]
Binary thinning:
\[ B_i\mid x_i \sim \mathrm{Bernoulli}(r_i/c), \qquad A=\sum_i B_i. \]
Useful fraction:
\[ \eta_{\mathrm{rej},N} = \frac{A}{N}, \qquad \mathbb E[\eta_{\mathrm{rej},N}\mid x_{1:N}] = \frac{1}{Nc}\sum_i r_i. \]
Large \(N\):
\[ \eta_{\mathrm{rej}} \to \frac{Z}{c}, \qquad Z=\mathbb E_q r. \]
Compare to IS:
\[ \eta_{\mathrm{IS}} = \frac{Z^2}{\mathbb E_q r^2}, \qquad \eta_{\mathrm{rej}} = \frac{Z}{c}. \]
von Neumann, “Various techniques used in connection with random digits”, NBS Appl. Math. Ser., 1951.
Function surrogate:
\[ q_\phi(x)\approx p(x). \]
Classical: Laplace/Gaussian, polynomial or spline quadrature.
Gaussian quadrature: exact for polynomial test functions.
Normalizing flow: invertible map of simple \(q_0\).
\[ z\sim q_0,\qquad x=T_\phi(z). \]
If approximate: reweight samples, \(\bar w_i\propto p(x_i)/q_\phi(x_i)\).
Problem: missed support stays invisible.
Quadrature:
\[ \int g(x)\rho(x)\,dx \approx \sum_i \lambda_i g(x_i). \]
Exact on a chosen function class.
Family:
\[ q_\phi\in\mathcal Q. \]
Flow:
\[ q_\phi=T_{\phi\#}q_0. \]
Fit:
\[ q_\phi = \arg\min_{q\in\mathcal Q} \mathrm{KL}(q\|p/Z). \]
\[ \mathrm{KL}(q\|p/Z) = \int q(x)\log\frac{q(x)}{p(x)/Z} \,dx. \]
Zero support:
\[ q_\phi(x)=0 \Rightarrow \text{no samples there}. \]
Laplace, “Memoire sur les approximations des formules qui sont fonctions de tres grands nombres”, Mem. Acad. R. Sci. Paris, 1785. Jordan et al., “An Introduction to Variational Methods for Graphical Models”, Mach. Learn., 1999. Rezende et al., “Variational Inference with Normalizing Flows”, ICML, 2015.
Global proposals might miss probability mass. Solution: move locally.
Define only a proposal kernel:
\[ x_{t+1}^\star\sim q(\cdot\mid x_t), \qquad q(x_{t+1}^\star\mid x_t). \]
Add Metropolis-Hastings (MH) correction (like rejection sampling):
\[ \alpha(x_t,x_{t+1}^\star) = \min\!\left( 1, \frac{p(x_{t+1}^\star)q(x_t\mid x_{t+1}^\star)} {p(x_t)q(x_{t+1}^\star\mid x_t)} \right). \]
\(\alpha\): acceptance probability.
Current state:
\[ x_t=x. \]
Proposal:
\[ y\sim q(\cdot\mid x). \]
Proposal ratio:
\[ r(x,y) = \frac{p(y)q(x\mid y)} {p(x)q(y\mid x)}. \]
Accept:
\[ \alpha(x,y)=\min(1,r(x,y)). \]
If not accepted, add the current value again:
\[ x_{t+1} = \begin{cases} y, & \text{with probability }\alpha(x,y),\\ x, & \text{otherwise}. \end{cases} \]
Transition:
\[ K(x,dy) = q(y\mid x)\alpha(x,y)dy +\rho(x)\delta_x(dy). \]
Sample measure:
\[ \widehat\Pi_N = \frac1N\sum_{t=1}^N\delta_{x_t}. \]
Metropolis et al., “Equation of State Calculations by Fast Computing Machines”, J. Chem. Phys., 1953. Hastings, “Monte Carlo Sampling Methods Using Markov Chains and Their Applications”, Biometrika, 1970.
Pairwise detailed balance: no net local flow.
\[ p(x)K(x,dy) = p(y)K(y,dx). \]
Integrate all balanced pairs: target density invariant.
\[ pK=p. \]
If the chain is ergodic, it converges to \(p\).
Then one long path gives posterior averages.
Target measure:
\[ \Pi(dx)=p(x)\,dx. \]
For any event \(A\):
\[ \begin{aligned} (\Pi K)(A) &= \int_\Omega p(x)K(x,A)\,dx\\ &= \int_\Omega\int_A p(x)K(x,dy)\,dx\\ &\overset{\mathrm{balance}}{=} \int_A\int_\Omega p(y)K(y,dx)\,dy\\ &= \int_A p(y)\,dy = \Pi(A). \end{aligned} \]
For the Markov chain:
\[ \begin{gathered} X_0=x_0,\qquad X_t\sim K^t(x_0,\cdot),\\ K^t(x_0,A) = \mathbb P(X_t\in A\mid X_0=x_0),\\ \text{balance} \Rightarrow \Pi K=\Pi \quad\text{(invariant target)},\\ \text{ergodicity conditions:}\\ \Pi\text{-irreducible},\quad \text{aperiodic},\quad \text{positive Harris}\\ \Rightarrow \lim_{t\to\infty}K^t(x_0,A)=\Pi(A) \quad\text{(independent of start)},\\ \Rightarrow \lim_{N\to\infty} \frac1N\sum_{t=1}^N f(X_t) = \int f\,d\Pi . \end{gathered} \]
Metropolis et al., “Equation of State Calculations by Fast Computing Machines”, J. Chem. Phys., 1953. Birkhoff, “Proof of the Ergodic Theorem”, Proc. Natl. Acad. Sci. U.S.A., 1931.
We already know: the transition kernel preserves the target.
\[ p K=p. \]
But local random-walk moves are not efficient:
\[ Y=X+\sigma\xi, \qquad \alpha(x,y)= \min\!\left\{1,\frac{p(y)}{p(x)}\right\}. \]
→ We want the proposal distribution to follow nearby probability mass.
Balance implies the invariant target:
\[ p(x)K(x,dy) = p(y)K(y,dx) \Rightarrow pK=p. \]
\[ Y=X+\sigma_d\xi, \qquad \sigma_d=\ell/\sqrt d. \]
For product targets:
\[ p_d(x)=\prod_{i=1}^d p_1(x_i), \qquad I=\mathbb E[(\partial_x\log p_1(X))^2]. \]
Then
\[ \bar\alpha_d(\ell) \to 2\Phi\!\left(-\frac{\ell\sqrt I}{2}\right), \qquad \ell_{\mathrm{opt}}: \bar\alpha\approx0.234. \]
So useful local moves require:
\[ \text{follow mass} \quad+\quad \text{small correction}. \]
Roberts et al., “Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms”, Ann. Appl. Probab., 1997.
Solution: gradient drift.
Drift the proposal toward higher density:
\[ \begin{aligned} Y &= X +\frac{h}{2}\nabla\log p(X) +\sqrt{h}\,\xi,\\ \xi&\sim N(0,I). \end{aligned} \]
After MH correction: MALA = Metropolis-adjusted Langevin algorithm.
Forward proposal:
\[ q_h(y\mid x) = N\!\left(y;\, x+\frac h2\nabla\log p(x),\, hI\right). \]
Balance compares flows:
\[ p(x)q_h(y\mid x) \quad\text{vs.}\quad p(y)q_h(x\mid y). \]
Therefore:
\[ \alpha(x,y) = \min\left( 1, \frac{p(y)q_h(x\mid y)} {p(x)q_h(y\mid x)} \right). \]
Not symmetric: contains drift \(\nabla\log p(x)\).
\[ q_h(y\mid x)\ne q_h(x\mid y) \]
Drift and diffusion flows cancel:
\[ \begin{gathered} dX_t = \frac12\nabla\log p(X_t)\,dt+dW_t,\\ \rho_t=\text{density}(X_t), \qquad W_t=\text{Brownian motion}. \end{gathered} \]
Probability current \(J_t\): net mass flux.
\[ J_t := \frac12\rho_t\nabla\log p - \frac12\nabla\rho_t, \qquad \partial_t\rho_t=-\nabla\cdot J_t. \]
At \(\rho_t=p\): \(J_t=\frac12p\nabla\log p-\frac12\nabla p=0\).
Rossky et al., “Brownian dynamics as smart Monte Carlo simulation”, J. Chem. Phys., 1978.
Remaining limitation: how far can one gradient be extrapolated?
Useful step length depends on how that gradient changes with \(x\).
Use \(U(x)=-\log p(x)\). Around the MAP point:
\[ x_{\mathrm{MAP}} := \underset{x}{\operatorname{arg\,max}}\,p(x) = \underset{x}{\operatorname{arg\,min}}\,U(x), \qquad \nabla U(x_{\mathrm{MAP}})=0. \]
The simplest local model: gradient is linear in \(\delta x\) with slope 1 (unit local scale).
\[ \begin{gathered} \delta x=x-x_{\mathrm{MAP}},\\ \nabla U(x_{\mathrm{MAP}}+\delta x) \approx \delta x. \end{gathered} \]
This is the Taylor / Gaussian local model.
Potential:
\[ U(x)=-\log(p(x)). \]
MAP point / minimum:
\[ x_{\mathrm{MAP}} = \underset{x}{\operatorname{arg\,max}}\,p(x) = \underset{x}{\operatorname{arg\,min}}\,U(x). \]
\[ \nabla U(x_{\mathrm{MAP}})=0, \qquad \delta x=x-x_{\mathrm{MAP}}. \]
Unit local scale:
\[ \nabla U(x_{\mathrm{MAP}}+\delta x) \approx \delta x. \]
Taylor at a minimum:
\[ U(x_{\mathrm{MAP}}+\delta x) \approx U(x_{\mathrm{MAP}}) + \frac12\|\delta x\|^2. \]
Local Gaussian:
\[ p_{\mathrm{loc}}(\delta x) \propto \exp\!\left(-\frac12\|\delta x\|^2\right). \]
If \(\nabla U\) changes like \(\delta x\), the measured gradient stays predictive over a local step.
Easy case for:
If the local model is anisotropic, “slopes” \(H\) differ in different directions:
\[ \nabla U(x_{\mathrm{MAP}}+\delta x)\approx H\delta x, \qquad H:=\nabla^2 U(x_{\mathrm{MAP}}). \]
with \(H\) the Hessian of \(U\), only diagonal if slopes are aligned with the axes.
Preconditioning problem with \(G = H^{-1}\) makes it look isotropic:
\[ G\approx H^{-1} \qquad\Rightarrow\qquad GH\approx I. \]
Equivalent view: change the distance scale.
Raw gradient change:
\[ \delta x\mapsto H\delta x. \]
Preconditioned change:
\[ \delta x\mapsto GH\delta x. \]
Nice units:
\[ GH\approx I. \]
Distance view:
same points, different unit lengths.
Choose a distance rule at each position:
\[ \|\delta x\|_{g(x)}^2 = (\delta x)^\top g(x)\,\delta x, \qquad g(x)\ \text{is symmetric positive definite (SPD)}. \]
The inverse \(G(x)=g(x)^{-1}\) scales gradients and noise.
If the distance rule changes with \(x\): Riemannian metric.
Infinitesimal distance:
\[ ds^2=dx^\top g(x)\,dx. \]
Inverse distance rule:
\[ G(x)=g(x)^{-1}. \]
Gradient in localized units:
\[ G(x)\nabla U(x). \]
Local quadratic prototype:
\[ g(x)\approx \nabla^2 U(x) \quad \text{if SPD.} \]
Otherwise use SPD proxies:
Fisher, Gauss–Newton, covariance, regularized Hessian.
Riemann, “Ueber die Hypothesen, welche der Geometrie zu Grunde liegen”, Abh. K. Ges. Wiss. Gottingen, 1868.