Suppose that \(X \sim F_X ( x )\) and let \(Y = g (X)\).
How do we find \(F_Y ( y )\) and \(f_Y ( y )\) ?
von Neumann
Given a uniform \(U\), how do we find \(X= g(U)\)?
In the bivariate case \((X,Y) \rightarrow (U,V)\).
We need to find \(f_{(U,V)} ( u , v )\) from \(f_{X,Y}(x,y)\)
Applications: Simulation, MCMC and PF.
Transformations
The cdf identity gives \[
F_Y ( y) = \mathbb{P} ( Y \leq y ) = \mathbb{P} ( g( X) \leq y )
\]
Hence if the function \(g ( \cdot )\) is monotone we can invert to get
\[
F_Y ( y ) = \int_{ g( x) \leq y } f_X ( x ) dx
\]
If \(g\) is increasing \(F_Y ( y ) = P( X \leq g^{-1} ( y ) ) = F_X ( g^{-1} ( y ) )\)
If \(g\) is decreasing \(F_Y ( y ) = P( X \geq g^{-1} ( y ) ) = 1 - F_X ( g^{-1} ( y ) )\)
Transformation Identity
Theorem 1: Let \(X\) have pdf \(f_X ( x)\) and let \(Y=g(X)\). Then if \(g\) is a monotone function we have
\[
f_Y ( y) = f_X ( g^{-1} ( y ) ) \left | \frac{ d}{dy} g^{-1} ( y ) \right |
\] There’s also a multivariate version of this that we’ll see later.
Suppose \(X\) is a continuous rv, what’s the pdf for \(Y = X^2\)?
Let \(X \sim N ( 0 ,1 )\), what’s the pdf for \(Y = X^2\)?
Probability Integral Transform
theorem Suppose that \(U \sim U[0,1]\), then for any continuous distribution function \(F\), the random variable \(X= F^{-1} (U)\) has distribution function \(F\).
Remember that for \(u \in [0,1]\), \(\mathbb{P} \left ( U \leq u \right ) = u\), so we have
\[
\mathbb{P} \left (X \leq x \right )= \mathbb{P} \left ( F^{-1} (U) \leq x \right )= \mathbb{P} \left ( U \leq F(x) \right )=F(x)
\] Hence, \(X = F_X^{-1}(U)\).
Normal
Sometimes thare are short-cut formulas to generate random draws
Normal \(N(0,I_2)\): \(x_1,x_2\) uniform on \([0,1]\) then \[
\begin{aligned}
y_1 = & \sqrt{-2\log x_1}\cos(2\pi x_2)\\
y_2 = & \sqrt{-2\log x_1}\sin(2\pi x_2)
\end{aligned}
\]
Simulation and Transformations
An important application is how to transform multiple random variables?
Suppose that we have random variables:
\[
( X , Y ) \sim f_{ X , Y} ( x , y )
\] A transformation of interest given by: \[
U = g ( X , Y )
\; \; {\rm and} \; \;
V = h ( X , Y )
\]
The problem is how to compute \(f_{ U , V } ( u , v )\) ? Jacobian
\[
J = \frac{ \partial ( x , y ) }{ \partial ( u , v ) } = \left | \begin{array}{cc}
\frac{ \partial x }{ \partial u} & \frac{ \partial x }{ \partial v} \\
\frac{ \partial y }{ \partial u} & \frac{ \partial y }{ \partial v}
\end{array} \right |
\]
Bivariate Change of Variable
Theorem: (change of variable)
\[
f_{ U , V } ( u , v ) = f_{ X , Y} ( h_1 ( u , v ) , h_2 ( u , v ) )
\left | \frac{ \partial ( x , y ) }{ \partial ( u , v ) } \right |
\] The last term is the Jacobian.
This can be calculated in two ways.
\[
\left | \frac{ \partial ( x , y ) }{ \partial ( u , v ) } \right | =
1 / \left | \frac{ \partial ( u , v ) }{ \partial ( x , y ) } \right |
\]
So we don’t always need the inverse transformation \(( x , y ) = ( g^{-1} ( u , v ) , h^{-1} ( u , v ) )\)
MCMC algorithms generate a Markov chain \[
\left\{ \theta ^{\left( g\right)}\right\} _{g=1}^{G}
\] whose stationary distribution is \(p\left( \theta \mid X,Y\right)\). Thus, the key to Bayesian inference is simulation rather than optimization.
Transition matrix \(P = \{p_{ij}\}_{i,j=1}^n\) is column stochastic, i.e. \(p_{ij} \ge 0\) and \(\sum_j p_{ij} = 1\).
Let \(\pi_k \in R^n\) be the distribution of the random variable \(X_k\) at time \(k\), \(\pi_{ki} = P(X_k=i)\)
\(\pi_{k+1} = P\pi_k\).
The limiting distribution is \(\pi = P\pi\), i.e. \(\pi\) is an eigenvector of \(P\) with eigenvalue 1.
The limiting distribution is unique and it is the first eigenvector of \(P\). The limiting distribution is also called a stationary distribution.
MCMC for Discrete Random Variables
Power method \(\pi^{k+1} = P^k\pi^0\) to \(\pi\) (Google’s PageRank).
Perron-Frobenius theorem: if \(p_{ij}>0\) then \(\pi\) always exists and it is unique and \[
\Vert \pi^k - \pi\Vert \le c |\lambda_2|^k,
\]
\(\lambda_2\) is the second largest eigenvalue of \(P\) (the first eigenvalue is always 1) - \(\pi\) as the average time of visiting vertex \(i\) under random walk.
The mixing time of the Markov chain is given by \[
T=\dfrac{1}{\log(1/\lambda_2)}
\] It is roughly, number of steps over which deviation from equilibrium distribution decreases by factor \(e\).
The reverse problem, is how to construct the transition probabilities so that the limiting distribution is \(\pi\).
For example, it is easy to see that any row-stochastic matrix which is symmetric has a uniform limiting distribution, i.e. \(\pi = e\), \(e_i = 1/n\), where \(n\) is the number of vertices.
\(e\) is an eigenvector of \(P^T\), since \(P^T\) is row-stochastic, thus sum of each row equals to 1, i.e. \(P^T e = e\).
General case balance condition: \[
\pi_i p_{ij} = \pi_j p_{ji}, ~ i,j=1,\ldots,n,
\] from which follows that \[
\sum_i p_{ij}\pi_i = \pi_j \sum_i p_{ji} = \pi_j,~j=1,\ldots,n.
\]
How to Construct Transition Probabilities?
Thus we need to find a Markov chain in which the transition probabilities satisfy the detailed balance condition (symmetry). Another way to write the condition is \[
\dfrac{p_{ij}}{p_{ji}} = \dfrac{\pi_j}{\pi_i}.
\] Say we have some “teaser” transition probabilities \(p_{ij}^0\) and we assume those are symmetric. We want to find new probabilities \(p_{ij} = p_{ij}^0b_{ij}, ~ i\ne j\) and \(p_{ii} = 1 - \sum_{j:~j\ne i} p_{ij}\). We need to choose \(b_{ij}\) is such a way so that \[
\dfrac{p_{ij}}{p_{ji}} = \dfrac{p_{ij}^0b_{ij}}{p_{ji}^0b_{ji}} = \dfrac{b_{ij}}{b_{ji}} = \dfrac{\pi_j}{\pi_i}
\] We can choose \[
b_{ij} = F\left(\dfrac{\pi_j}{\pi_i}\right)
\]
For the detailed balance condition to be satisfied, we need to choose \(F: R_+ \rightarrow [0,1]\) such that \[
\dfrac{F(z)}{F(1/z)} = z.
\]
How to Construct Transition Probabilities?
An example of such a function is \(F(z) = \min(z,1)\). This leads to the Metropolis algorithm
When at state \(i\), draw from \(p^{0}_{ij}\) to generate next state \(j\)
Move to \(j\) with probability \(a_{ij}\), stay at \(i\) otherwise
When move happens, we say step is accepted, otherwise it is rejected. In a more general case (Metropolis-Hastings), when \(p^0\) is not symmetric, we calculate \[
a_{ij} = \min\left(1,\dfrac{\pi_i p^0_{ji}}{\pi_j p^0_{ji}}\right).
\]
Continious Case
In the continuous case we use variables \(S\) and \(T\) instead induces \(i\) and \(j\). The detailed balance condition becomes \[
\pi(S)p(T \mid S) = \pi(T)p(S \mid T).
\] Then \[
\int p(T \mid S)\pi(S)dS = \int p(S \mid T)\pi(T)dS = \pi(T)\int p(S\mid T)dS = \pi(T).
\]
The following conditions \[
\forall S, \forall T:\pi(T)\ne 0:p(T \mid S) >0
\] are sufficient to guarantee uniqueness of \(\pi(T)\).
Metropolis-Hastings Algorithm
Specifically, the MH algorithm repeats the following two steps \(G\) times: given \(\theta ^{\left( g\right) }\)\[
\begin{aligned}
&\text{Step 1. Draw }\theta^{\prime} \text{ from a proposal distribution,} p(\theta^{\prime}|\theta ^{(g)}) \\
&\text{Step 2. Accept } \theta^{\prime} \text{ with probability } \alpha \left(\theta ^{(g)},\theta^{\prime}\right) \text{,}
\end{aligned}
\] where \[
\alpha \left( \theta ^{(g)},\theta^{\prime}\right) =\min \left( \frac{\pi(\theta^{\prime})}{\pi (\theta ^{(g)})}\frac{q(\theta^{(g)}|\theta^{\prime})}{q(\theta^{\prime}|\theta ^{(g)})},1\right) \text{.}
\]
Metropolis-Hastings Algorithm
Implementation of the accept-reject step:
Draw a uniform random variable, \(U\sim U \left[ 0,1\right]\), and set \(\theta ^{\left( g+1\right) }=\theta^{\prime}\)
It is important to note that the denominator in the acceptance probability cannot be zero, provided the algorithm is started from a \(\pi\) - positive point since \(q\) is always positive. The MH algorithmonly requires that \(\pi\) can be evaluated up to proportionality.
The output of the algorithm, \(\left\{ \theta ^{\left( g\right) }\right\}_{g=1}^{\infty }\), is clearly a Markov chain. The key theoretical property is that the Markov chain, under mild regularity, has \(\pi \left( \theta\right)\) as its limiting distribution. We discuss two important specialcases that depend on the choice of \(q\).
Independence MH
One special case draws a candidate independently of the previous state, \(q(\theta^{\prime}|\theta ^{(g)})=q(\theta^{\prime})\). In this independence MH algorithm, the acceptance criterion simplifies to \[
\alpha \left( \theta ^{(g)},\theta^{\prime}\right) =\min \left( \frac{\pi (\theta^{\prime})}{\pi (\theta ^{(g)})}\frac{q(\theta ^{(g)})}{q(\theta ^{\prime})},1\right).
\] Even though \(\theta\) is drawn independently of the previous state, the sequence generated is not being independent, since \(\alpha\) depends on previous draws. The criterion implies a new draw is always accepted if target density ratio \(\pi (\theta^{\prime})/\pi (\theta ^{(g)})\), increases more than the proposal ratio, \(q(\theta ^{(g)})/q(\theta^{\prime})\). When this is not satisfied, a balanced coin is flipped to decide whether or not toaccept the proposal.
Random-walk Metropolis
Random-walk (RW) Metropolis is the polar opposite of the independence MH algorithm. It draws a candidate from the following RW model, \[
\theta^{\prime}=\theta ^{\left( g\right) }+\sigma \varepsilon _{g+1},
\] where \(\varepsilon _{t}\) is an independent, mean zero, and symmetric error term, typically taken to be a normal or \(t-\)distribution, and \(\sigma\) is a scaling factor. The algorithm must be tuned via the choice of \(\sigma\), the scaling factor. Symmetry implies that \[
q\left( \theta^{\prime}|\theta ^{\left( g\right) }\right) =q\left( \theta ^{\left( g\right) }|\theta^{\prime}\right),
\] and \[
\alpha \left( \theta ^{(g)},\theta^{\prime}\right) =\min \left( \pi (\theta^{\prime})/\pi (\theta ^{(g)}),1\right).
\]
Implementation
set.seed(7)mh =function(target, proposal, n, x0) { x = x0; p =length(x0) samples =matrix(NA, nrow=n, ncol=p) accept =rep(0, n)for (i in1:n) { x_new =proposal(x) a =min(1,target(x_new) /target(x))# print(c(x, x_new, a))if (runif(1) < a) {accept[i]=1; x = x_new} samples[i,] = x# print(accept[i]) }list(samples = samples, accept = accept)}
Apply
We apply MCMC to the weighted sum of two normal distributions. This sort of distribution is fairly straightforward to sample from, but let’s draw samples with MCMC.
Clifford-Hammersley (CH) theorem: high-dimensional joint distribution, \(p\left( \theta\mid X,Y\right)\), is completely characterized by a larger number of lower dimensional conditional distributions. Given this characterization, MCMC methods iteratively sample from these conditional distributions using standard sampling methods.
Example: bivariate posterior distribution \(p\left( \theta_1,\theta_{2} \mid X,Y \right)\), e.g. \(\theta_1\) as traditional static parameters and \(\theta_2\) as latent variables.
CH: \(p\left( \theta_1,\theta_2\right)\) is determined by \(p\left( \theta_{1} \mid \theta_{2}\right)\) and \(p\left( \theta_{2} \mid \theta_{1}\right)\), given \(p\left( \theta_{1},\theta_2\right)\), \(p\left( \theta_{1}\right)\) and \(p\left( \theta_{2}\right)\) have positive mass for all points
Besag formula
For any pairs \((\theta_{1},\theta_{2})\) and \((\theta_1^{\prime},\theta_2^{\prime})\), \[
\frac{p(\theta_1,\theta_2)}{p(\theta_1^{\prime},\theta_2^{\prime})}=\frac{p(\theta_1 \mid \theta_2^{\prime})p(\theta_2 \mid \theta_1)}{p\left( \theta_1^{\prime} \mid \theta_2^{\prime}\right) p\left( \theta_2^{\prime} \mid \theta_2\right) }.
\] The proof uses the fact that \(p\left( \theta_1,\theta_2\right)=p\left( \theta_2 \mid \theta_1\right) p\left( \theta_1\right)\) (applied to both \((\theta_1,\theta_2)\) and \((\theta_1^{\prime},\theta_2^{\prime})\)) and the Bayes rule: \[
p\left( \theta_1\right) =\frac{p\left( \theta_1 \mid \theta_2^{\prime }\right) p\left( \theta_2^{\prime}\right) }{p\left( \theta_2^{\prime} \mid \theta_1\right) }\text{.}
\]
Clifford-Hammersley: The general version
Partitioning a vector as \(\theta =\left( \theta_{1},\theta_{2},\theta_{3},\ldots ,\theta_{K}\right)\), then the general CH theorem states that \[
p\left( \theta _{i} \mid \theta _{-i}\right) := p\left( \theta _{i} \mid \theta_1,\theta_1,\ldots,\theta_{i-1},\theta _{i+1},...,\theta _{K}\right) ,
\] for \(i=1,...,K\), completely characterize \(p\left( \theta_1,\ldots,\theta_{K}\right)\).
Latent Variable Models
Fixed parameters \(\theta\), and latent variables, \(x\). In this case, then CH \[
p\left( \theta,x \mid y\right)
\] is completely characterized by \(p\left( \theta \mid x,y\right)\) and \(p\left( x \mid \theta ,y\right)\).
The distribution \(p\left( \theta \mid x,y\right)\) is the posterior distribution of the parameters, conditional on the observed data and the latent variables. Similarly, \(p\left( x \mid \theta,y\right)\) is the smoothing distribution of the latent variables.
Gibbs Sampler
The Gibbs sampler simulates multi-dimensional posterior distributions by iteratively sampling from the lower-dimensional conditional posteriors.
Unlike the previous MH algorithms, the Gibbs sampler updates the chain one component at a time, instead of updating the entire vector.
This requires either that the conditional posteriors distributions is discrete, or a recognizable distribution (e.g. Normal) for which standard sampling algorithms apply, or that resampling methods, such as accept-reject, can be used.
Gibbs Sampler
In the case of \(p\left( \theta_1,\theta_2\right)\), given current draws,\(\left( \theta_1^{\left( g\right) },\theta_2^{\left( g\right)}\right)\), the Gibbs sampler consists of \[
\begin{aligned}
\text{1. Draw }\theta_1^{\left( g+1\right) } &\sim &p\left( \theta _1|\theta_2^{\left( g\right) }\right) \\
\text{2. Draw }\theta_2^{\left( g+1\right) } &\sim &p\left( \theta_2|\theta_1^{\left( g+1\right) }\right) ,
\end{aligned}
\] repeating \(G\) times.
Gibbs Sampler
Transition kernel from state \(\theta\) to state \(\theta ^{\prime}\) is \[
p\left( \theta ,\theta^{\prime}\right) =p\left( \theta_{1}^{\prime}|\theta_2\right) p\left( \theta_2^{\prime}|\theta_{1}^{\prime}\right)
\]
The limiting (statiobnary) distribution \(\pi\), satisfies \[
\pi \left( \theta^{\prime}\right) =\int p\left( \theta ,\theta^{\prime}\right) d\theta.
\]
Gibbs Sampler
It is easy to verify that the stationary distribution of the Markov chain generated by the Gibbs sampler is the posterior distribution, \(\pi \left(\theta \right) =p\left( \theta_1,\theta_2\right)\): \[
\begin{aligned}
\int p\left( \theta ,\theta^{\prime}\right) d\theta &=&p\left( \theta_2^{\prime}|\theta_1^{\prime }\right) \int_{\theta_2}\int_{\theta_1}p\left( \theta_1^{\prime}|\theta_2\right) p\left( \theta_1,\theta_2\right) d\theta_{1}d\theta_2 \\
&=&p\left( \theta_2^{\prime}|\theta_1^{\prime}\right) \int_{\theta_2}p\left( \theta_1^{\prime}|\theta_2\right) p\left( \theta_{2}\right) d\theta_2 \\&=&p\left( \theta_2^{\prime}|\theta_1^{\prime}\right) p\left( \theta_{1}^{\prime}\right) =p\left( \theta_1^{\prime},\theta_2^{\prime}\right) =\pi \left( \theta^{\prime}\right) \text{.}
\end{aligned}
\]
Gibbs Sampler
The ergodic theorem holds: for a sufficiently integrable function \(l\) and for all starting points \(\theta ,\)\[
\underset{G\rightarrow \infty }{\lim }\frac{1}{G}\sum_{g=1}^{G}f\left(\theta ^{\left( g\right) }\right) =\int f\left( \theta \right) \pi \left(\theta \right) d\theta =E\left[ f\left( \theta \right) \right]
\] - Run for an initial length, often called the burn-in - Then a secondary sample of size \(G\) is created for Monte Carlo inference
Gibbs sample for Normal-Gamma
Code
# summary statistics of samplen <-30ybar <-15s2 <-3# sample from the joint posterior (mu, tau | data)mu <- tau <-rep(NA, 11000)T <-1000# burnintau[1] <-1# initialisationfor(i in2:11000){ mu[i] <-rnorm(n =1, mean = ybar, sd =sqrt(1/ (n * tau[i -1]))) tau[i] <-rgamma(n =1, shape = n /2, scale =2/ ((n -1) * s2 + n * (mu[i] - ybar)^2))}mu <- mu[-(1:T)] # remove burnintau <- tau[-(1:T)] # remove burninhist(mu)
Code
hist(tau)
Hybrid chains
Given a partition of the vector \(\theta\) via CH, a hybrid MCMC algorithm updates the chain one subset at a time, either by direct draws (Gibbs steps) or via MH. - Thus, a hybrid algorithm combines the features of the MH algorithm and the Gibbs sampler, providing significant flexibility in designing MCMC algorithms for different models.
\(p\left( \theta_2|\theta_1\right)\) is recognizable and can be directly sampled.
\(p\left( \theta_1|\theta_2\right)\) can only be evaluated and not directly sampled.
The general hybrid algorithm is as follows. Given \(\theta_1^{\left(g\right) }\) and \(\theta_2^{\left( g\right) }\), for \(g=1,\ldots, G\), \[
\begin{aligned}
1.\text{ Draw }\theta_1^{\left( g+1\right) } &\sim &MH\left[ q\left( \theta_1|\theta_1^{\left( g\right) },\theta_2^{\left( g\right) }\right) \right] \\
2.\text{ Draw }\theta_2^{\left( g+1\right) } &\sim &p\left( \theta_2|\theta_1^{\left( g+1\right) }\right) .
\end{aligned}
\] In higher dimensional cases, a hybrid algorithm consists of any combination of Gibbs and Metropolis steps. Hybrid algorithms significantly increase the applicability of MCMC methods, as the only requirement is that the model generates posterior conditionals that can either be sampled or evaluated.
Hamiltonian Monte-Carlo
Add a momentum to the Metropolis Sampling.
Gradient of the un-normalized posterior to better navigate the surface defined by the posterior.
Momentum variable \(r\) for each model variable \(\theta\). \[
p(\theta,r) \propto \exp\left(L(\theta) - 0.5 r^Tr\right),
\]
\(L\) is the logarithm of the joint density of the variables of interest (negative potential energy)
\(0.5 r^Tr\) is the kinetic energy of the particle
\(p(\theta,r)\) is the total negative energy
Hamiltonian Monte-Carlo
We can simulate the evolution over time of the Hamiltonian dynamics of this system via the “leapfrog” integrator, which proceeds according to the updates \[
\begin{aligned}
r^{+/2} & = r + (\epsilon/2) \nabla_{\theta}L(\theta) \\
\theta^{+} & = \theta + (\epsilon/2) r^{+/2}\\
r^{+} & = r^{+/2} + (\epsilon/2)\nabla_{\theta}L(\theta^{+}).
\end{aligned}
\]
Hierarchical Bayesian Model
We consider the following model: \[\begin{align*}
\tau &\sim \mathrm{Gamma}(0.5, 0.5) \\
\lambda_d &\sim \mathrm{Gamma}(0.5, 0.5) \\
\beta_d &\sim \mathcal{N}(0, 20) \\
y_n &\sim \mathrm{Bernoulli}(\sigma((\tau \lambda \odot \beta)^T x_n))),
\end{align*}\]
\(\tau\) is a scalar global coefficient scale
\(\lambda\) is a vector of local scales
\(\beta\) is the vector of unscaled coefficients,
Code
import tensorflow_probability.substrates.jax as tfpimport numpy as npimport jax.numpy as jnptfd = tfp.distributionsimport jaximport matplotlib.pyplot as plt
Apply to Iris Data
Code
import pandas as pdimport matplotlib.pyplot as pltimport numpy as np iris = pd.read_csv("data/iris.csv")print(iris.shape)
Taking the logarithm of both sides, we get \[
\log \pi(z) = \log \pi(\theta) + \log \left| \frac{\partial T }{\partial z}(z) \right|
\] Use \(T(z)=e^z\), and \(\log|\frac{\partial T}{\partial z}(z)| = z\)
Change of Variables in Code
def unconstrained_joint_log_prob(x, y, theta): ndims = x.shape[-1] unc_tau, unc_lamb, beta = jnp.split(theta, [1, 1+ ndims]) unc_tau = unc_tau.reshape([]) # Make unc_tau a scalar tau = jnp.exp(unc_tau) ldj = unc_tau lamb = jnp.exp(unc_lamb) ldj += unc_lamb.sum()return joint_log_prob(x, y, tau, lamb, beta) + ldj
Let’s check out function
from functools import partialtarget_log_prob = partial(unconstrained_joint_log_prob, x, y)theta = np.r_[tau,lamb,beta]target_log_prob(theta)
y = iris$Species=="setosa"glm(y~iris$Sepal.Length, family=binomial)
Call: glm(formula = y ~ iris$Sepal.Length, family = binomial)
Coefficients:
(Intercept) iris$Sepal.Length
27.83 -5.18
Degrees of Freedom: 149 Total (i.e. Null); 148 Residual
Null Deviance: 191
Residual Deviance: 72 AIC: 76
Hierarchical Model
We will use the Normal-Normal model with kown variance and unknowm mean \[
\begin{aligned}
\mu \sim & \mathrm{Normal}(0,1) \\
x\mid \mu \sim & \mathrm{Normal}(\mu,1)
\end{aligned}
\]
MH is very simple and quite general. At the same time
too large step size \(\sigma\) leads to a large fraction of rejected samples, while too small \(\sigma\) makes very small steps, thus it takes long time to ‘explore the distribution’ (check this in demonstration!)
in high-dimensional spaces (very important use-case), MH explores the space very inefficiently because of it’s random-walk behavior. Guessing good direction in 1000 dimensions is incomparably harder than doing this for 2 dimensions!
MH can’t travel long distances (significantly larger than \(\sigma\)) between isolated local minimums
Sampling high-dimensional distributions with MH becomes very inefficient in practice. A more efficient scheme is called Hamiltonian Monte Carlo (HMC).
