Bayesian modelling

Monte Carlo methods and Markov chains

Léo Belzile

Last compiled Tuesday Sep 22, 2026

Bayesian inference beyond conjugate models

How to circumvent the problem of intractable posteriors?

  • simulation-based methods: accept-reject, Markov chain Monte Carlo, particle filters, etc.
  • deterministic methods: (integrated nested) Laplace approximations, variational Bayes, expectation propagation, etc.

We focus on Monte Carlo methods.

Simulation algorithms: inversion method

If \(F\) is an absolutely continuous distribution function, then \[F(X) \sim \mathsf{unif}(0,1).\] The inversion method consists in applying the quantile function \(F^{-1}\) to \(U \sim \mathsf{unif}(0,1)\), viz. \[F^{-1}(U) \sim X.\]

Inversion method for truncated distributions

Consider a random variable \(Y\) with distribution function \(F\).

If \(X\) follows the same distribution as \(Y\), but restricted over the interval \([a,b]\), then \[\Pr(X \leq x) = \frac{F(x) - F(a)}{F(b)-F(a)}, \qquad a \leq x \leq b,\]

Therefore, \[F^{-1}[F(a) + \{F(b)-F(a)\}U] \sim X.\]

Fundamental theorem of simulation

Consider a \(d\)-variate random vector \(\boldsymbol{X},\) independent of \(U \sim \mathsf{unif}(0,1)\) and \(c>0\). If \((\boldsymbol{X}, U)\) is uniformly distributed on the set \[\begin{align*} \mathcal{A}_{f}=\{(\boldsymbol{x}, u): 0 \leq u \leq c f(\boldsymbol{x})\}, \end{align*}\] then \(\boldsymbol{X}\) has density \(f(\boldsymbol{x}).\)

  • \(f\) is the marginal density of \(\boldsymbol{X}\) since \(f(\boldsymbol{x}) = \int_0^{f(\boldsymbol{x})} \mathrm{d} u.\)
  • If we can simulate uniformly from \(\mathcal{A}_{f},\) then, we can discard the auxiliary variable \(u.\) See Devroye (1986), Theorem 3.1.

Fundamental theorem of simulation in picture

Figure 1: Illustration of the fundamental theorem of simulation. All points in blue below the density curve belong to \(\mathcal{A}_f.\)

Simulation algorithms: accept-reject

  • Target: sample from density \(p(x)\) (hard to sample from)
  • Proposal: find a density \(q(x)\) with nested support, \(\mathrm{supp}(p) \subseteq \mathrm{supp}(q)\), such that for all \(x \in \mathrm{supp}(p)\), \[\frac{p(x)}{q(x)} \leq C, \quad C \geq 1.\] Uses the fundamental theorem of simulation by finding an envelope \(Cq(x)\) that is over the density \(p(x).\)

Rejection sampling algorithm

  1. Generate \(X\) from proposal with density \(q(x)\).
  2. Compute the ratio \(R \gets p(X)/ q(X)\).
  3. If \(CU \leq R\) for \(U \sim \mathsf{unif}(0,1)\), return \(X\), else go back to step 1.

Accept-reject illustration

Figure 2: Target density (full) and scaled proposal density (dashed): the vertical segment at \(x=1\) shows the percentage of acceptance for a uniform slice under the scaled proposal, giving an acceptance ratio of 0.58.

Remarks on rejection sampling

  • Acceptance rate is \(1/C\):
    • we need on average \(C\) draws from \(q\) to get one from \(p\).
  • \(q\) must be more heavy-tailed than \(p\)
    • e.g., \(q(x)\) Student-\(t\) for \(p(x)\) Gaussian.
  • \(q\) should be cheap and easy to sample from!

Designing a good proposal density

Good choices must satisfy the following constraints:

  • pick a family \(q(x)\) so that \[C = \mathrm{argmax}_x{p(x)}/{q(x)}\] is as close to 1 as possible.
  • you can use numerical optimization with target \[\log R(x) =\log p(x) - \log q(x)\] to find the mode \(x^\star\) and the upper bound \(C = R(x^\star).\)
  • if \(q(x; \boldsymbol{\theta})\) is a function of the parameter \(\boldsymbol{\theta}\), search for \[\mathrm{argmin}_{\boldsymbol{\theta}} \mathrm{argmax}_x \{\log p(x) - \log q(x; \boldsymbol{\theta})\}\]

Truncated Gaussian via accept-reject

Consider sampling \(Y \sim \mathsf{Gauss}(\mu, \sigma^2)\), but truncated in the interval \((a, b)\). The target density is \[\begin{align*} p(x; \mu, \sigma, a, b) = \frac{1}{\sigma}\frac{\phi\left(\frac{x-\mu}{\sigma}\right)}{\Phi(\beta)-\Phi(\alpha)}. \end{align*}\] for \(\alpha= (a-\mu)/\sigma\) and \(\beta = (b-\mu)/\sigma\). where \(\phi(\cdot), \Phi(\cdot)\) are respectively the density and distribution function of the standard Gaussian distribution.

Accept-reject (crude version)

  1. Simulate \(X \sim \mathsf{Gauss}(\mu, \sigma^2)\)
  2. reject any draw if \(X < a\) or \(X> b\).

The acceptance rate is \(C^{-1} = \{\Phi(\beta) - \Phi(\alpha)\}\)

# Standard Gaussian truncated on [0,1]
candidate <- rnorm(1e5)
trunc_samp <- candidate[candidate >= 0 & candidate <= 1]
# Acceptance rate
length(trunc_samp)/1e5
[1] 0.33972
# Theoretical acceptance rate
pnorm(1)-pnorm(0)
[1] 0.3413447

Accept-reject for truncated Gaussian

Since the Gaussian is a location scale family, the inversion method gives \[\begin{align*} X \sim \mu + \sigma\Phi^{-1}\left[\Phi(\alpha) + \{\Phi(\beta)-\Phi(\alpha)\}U\right] \end{align*}\]

We however need to evaluate \(\Phi\) numerically (no closed-form expression).

The method fails for rare event simulation because the computer returns

  • \(\Phi(x) = 0\) for \(x \leq -39\)
  • \(\Phi(x)=1\) for \(x \geq 8.3\),

implying that \(a \leq 8.3\) for this approach to work (Botev & L’Écuyer, 2017).

Simulating tails of Gaussian variables

We consider simulation from a standard Gaussian truncated above \(a>0\)

Write the density of the truncated Gaussian as (Devroye, 1986, p. 381)\[f(x) = \frac{\exp(-x^2/2)}{\int_{a}^{\infty}\exp(-z^2/2)\mathrm{d} z} =\frac{\exp(-x^2/2)}{c_1}.\]

Note that, for \(x \geq a\), \(x/a > 1\) and \[c_1f(x) \leq \frac{x}{a}\exp\left(-\frac{x^2}{2}\right)= a^{-1}\exp\left(-\frac{a^2}{2}\right)g(x);\] where \(g(x)=x\exp\{(a^2-x^2)/2\}\mathsf{1}_{x \ge a}\) is the density of a Rayleigh variable shifted by \(a\).

The constant \(C= \exp(-a^2/2)(c_1a)^{-1}\) approaches 1 quickly as \(a \to \infty\) (asymptotically optimal).

Accept-reject: truncated Gaussian with Rayleigh

The shifted Rayleigh has distribution function \[G(x) = 1-\exp\{(a^2-x^2)/2\}, x \geq a.\]

Marsaglia algorithm

  1. Generate a shifted Rayleigh above \(a\), \(X \gets \{a^2 - 2\log(U)\}^{1/2}\) for \(U \sim \mathsf{unif}(0,1)\)
  2. Accept \(X\) if \(XV \leq a\), where \(V \sim \mathsf{unif}(0,1)\).

For sampling on \([a,b]\) with \(a\) very large, propose from a Rayleigh truncated above at \(b\) (Botev & L’Écuyer, 2017).

a <- 8.3
niter <- 1000L
X <- sqrt(a^2 + 2*rexp(niter))
samp <- X[runif(niter)*X <= a]

Accept if \(U \leq \frac{f(x)}{g(x)} \times \frac{1}{C} = \frac{\exp(-x^2/2)/c_1}{x\exp\{(a^2-x^2)/2\}} \times \exp(a^2/2)c_1a = \frac{a}{x}.\)

Monte Carlo integration

Calculate \(\mathsf{E}_{X}\{g(X)\}\) of \(X\) with density function \(p\) on \(\mathcal{X}.\)

Monte Carlo integration proceeds by drawing \(B\) independent samples \(x_1, \ldots, x_B\) from \(p\) on \(\mathcal{X}\) and evaluating the empirical average of \(g,\) \[\begin{align*} \mathsf{E}\{g(X)\} = \int_{\mathcal{X}} g(x) p(x) \mathrm{d} x \approx \widehat{\mathsf{E}}\{g(X)\}=\frac{1}{B}\sum_{b=1}^B g(x_b). \end{align*}\]

Importance sampling

We can calculate the integral \(\mathsf{E}\{g(X)\} = \int_{\mathcal{X}} g(x) p(x) \mathrm{d}x\) by considering instead draws from a density \(q(\cdot)\) if \(\mathcal{X} \subseteq \mathrm{supp}(q).\) Then, \[\begin{align*} \mathsf{E}\{g(X)\} = \int_{\mathcal{X}} g(x) \frac{p(x)}{q(x)} q(x) \mathrm{d} x. \end{align*}\]

Monte Carlo estimator for importance sampling

An alternative Monte Carlo estimator, which is biased but has lower variance, is obtained by drawing independent \(x_1, \ldots, x_B\) from \(q\) and taking instead the weighted average of \[\begin{align*} \widetilde{\mathsf{E}}\{g(X)\} =\frac{B^{-1} \sum_{b=1}^B w_b g(x_b) }{B^{-1}\sum_{b=1}^B w_b}. \end{align*}\] with weights \(w_b = p(x_b)/q(x_b).\) The latter equal 1 on average, so one could omit the denominator without harm.

Importance sampling for the variance of a beta distribution

Consider \(X \sim \mathsf{beta}(\alpha, \alpha)\) for \(\alpha > 1\), with expectation \(\mathsf{E}(X)=0.5\) since the density is symmetric.

The variance is \[\mathsf{Va}(X) = \mathsf{E}\{(X - 0.5)^2\}=\frac{1}{4 \cdot (2\alpha+1)}.\]

We consider an equiweighted (bimodal) mixture of \(\mathsf{beta}(\alpha, 3\alpha)\) and \(\mathsf{beta}(3\alpha, \alpha)\) as importance sampling density.

Importance sampling for variance of a beta distribution

Code for importance sampling of beta variance

B <- 2e4L; alpha <- 1.5; factor <- 3
# Mode at the mean 0.5
X0 <- rbeta(n = B, alpha, alpha)
px <- function(x){dbeta(x, alpha, alpha)}
# Importance sampling density - mixture of two betas (alpha, factor*alpha)
X1 <- ifelse(
  runif(B) < 0.5,
  rbeta(B, alpha, factor*alpha), 
  rbeta(B, factor*alpha, alpha))
qx <- function(x){
  0.5*dbeta(x, alpha, factor*alpha) + 0.5*dbeta(x, factor*alpha, alpha)}
# Function to integrate - gives variance of a symmetric beta distribution
g <- function(x){(x - 0.5)^2}
# Weights for importance sampling
w <- px(X1)/qx(X1)

Code for importance sampling of beta variance (continued)

# Monte Carlo integration
mc_est <- mean(g(X0))
mc_var <- var(g(X0))/B
# Importance sampling weighted mean and variance
is_est <- weighted.mean(g(X1), w = w) # equivalent to mean(g(X1)*w)/mean(w)
is_var <- sum(w^2*(g(X1) - is_est)^2)/ (sum(w)^2)
# True value for the beta variance
th_est <- 1/(4*(2*alpha+1))
# Point estimates and differences
round(c(true = th_est, "monte carlo" = mc_est, "importance sampling" = is_est),4)
               true         monte carlo importance sampling 
             0.0625              0.0627              0.0618 
# Ratio of std. errors for means
mc_var/is_var # value > 1 means that IS is more efficient
[1] 1.119518

Markov chains

Plain ordinary Monte Carlo is great, but few algorithms are generic enough to be useful in complex high-dimensional problems.

We will instead typically build Markov chains that target an invariant stationary distribution.

Caveats?

Markov chain Monte Carlo methods generate correlated draws.

Questions:

  1. can we use them as ordinary independent samples?
  2. what is the price to pay?

We need to do a little theoretical detour to answer these questions.

Stationarity and Markov property

A stochastic process is

  • (strongly) stationary if the distribution of \(\{X_1, \ldots, X_t\}\) is the same as that of \(\{X_{n+1}, \ldots X_{t+n}\}\) for any value of \(n\) and given \(t\).
  • weakly stationary if \(\mathsf{E}(X_t) = \mu\) for all \(t\), and \(\mathsf{Cov}(X_t, X_{t+h}) = \gamma_h\) does not depend on \(t\).
  • Markov if it satisfies the Markov property: given the current state of the chain, the future only depends on the current state and not on the past.

Strong stationarity implies weak stationarity.

Law of large number (ergodic theorem)

Let \(\{Y_t\}\) is a weakly stationary sequence with mean \(\mathsf{E}(Y_t)=\mu\) and \(\gamma_h = \mathsf{Cov}(Y_t, Y_{t+h})\). Then, if the autocovariance series is convergent, meaning \[\sum_{t=0}^\infty |\gamma_h| < \infty,\] then \(\{Y_t\}\) is ergodic for the mean and \(\overline{Y} \stackrel{\mathrm{p}}{\to} \mu\).

Ergodicity means that two segments of a time series far enough apart act as independent.

Ergodicity and transformations

The ergodic theorem is a law of large numbers for stochastic processes that allows for serial dependence between observations, provided the latter is not too large.

Any transformation \(g(\cdot)\) of a stationary and ergodic process \(\{Y_t\}\) retains the properties, so \(\overline{g} = T^{-1} \sum_{t=1}^T g(Y_t) \to \mathsf{E}\{g(Y_t)\}\) as \(T \to \infty.\)

Autoregressive process of order 1

Consider a first-order autoregressive process, or \(\mathsf{AR}(1)\),

\[Y_t = \mu + \phi(Y_{t-1} - \mu) + \varepsilon_t,\] where

  • \(\phi\) is the lag-one correlation,
  • \(\mu\) the global mean
  • \(\varepsilon_t\) is an iid innovation with mean zero and variance \(\sigma^2\)

If \(|\phi| < 1\), the process is stationary, otherwise the variance increases with \(t\).

Unconditional moments via tower law

If the process is weakly stationary, then \(\mathsf{E}_{Y_{t}}(Y_t)=\mathsf{E}_{Y_{t-1}}(Y_{t-1})\) \[\begin{align*} \mathsf{E}_{Y_{t}}(Y_t) &= \mathsf{E}_{Y_{t-1}}\left\{\mathsf{E}_{Y_{t} \mid Y_{t-1}}(Y_t)\right\} \\&= \mu(1-\phi) + \phi\mathsf{E}_{Y_{t-1}}(Y_{t-1}) \end{align*}\] and so the unconditional mean is \(\mu\). For the variance, we have \[\begin{align*} \mathsf{Va}_{Y_{t}}(Y_t) &= \mathsf{E}_{Y_{t-1}}\left\{\mathsf{Va}_{Y_{t} \mid Y_{t-1}}(Y_t)\right\} + \mathsf{Va}_{Y_{t-1}}\left\{\mathsf{E}_{Y_{t} \mid Y_{t-1}}(Y_t)\right\}\\ & = \sigma^2 + \mathsf{Va}_{Y_{t-1}}\left\{\mu + \phi(Y_{t-1} - \mu)\right\} \\&= \sigma^2 + \phi^2 \mathsf{Va}_{Y_{t-1}}(Y_{t-1}). \end{align*}\] and the unconditional variance is \[\mathsf{Va}(Y_t) = \tau^2 = \sigma^2/(1-\phi^2).\]

Autocovariance and Markov property

The covariance at lag \(h\), in terms of innovations, gives \[\begin{align*} \gamma_h = \mathsf{Co}(Y_t, Y_{t-h}) = \mathsf{Va}(\phi Y_{t-1}, Y_{t-h}) + \mathsf{Va}(\varepsilon_t, Y_{t-h}) = \phi \gamma_{h-1} \end{align*}\] since \(\varepsilon_t\) is independent of the past. We thus find \[\rho_h = \phi^h, \qquad \gamma_h = \phi^h\mathsf{Va}(Y_t).\]

The \(\mathsf{AR}(1)\) process is first-order Markov since the conditional distribution \(p(Y_t \mid Y_{t-1}, \ldots, Y_{t-p})\) equals \(p(Y_t \mid Y_{t-1}).\)

Variance of sample average

Intuitively, a sample of correlated observations carries less information than an independent sample of draws.

The variance of the sample average is \[\begin{align*} \mathsf{Va}\left(\overline{Y}_T\right) &= \frac{1}{T^2}\sum_{t=1}^T \sum_{s=1}^T \mathsf{Co}(Y_t, Y_s) \\&= \frac{1}{T^2}\sum_{t=1}^T \mathsf{Va}(Y_t) + \frac{2}{T^2} \sum_{t=1}^{T-1}\sum_{s = t+1}^T \mathsf{Co}(Y_t, Y_s). \end{align*}\] Under independence and assuming stationarity, we get \(\mathsf{Va}\left(\overline{Y}_T\right)=\sigma^2/T.\)

Variance of sample average, redux

If the second moments are finite, the scaled limiting variance of the sample mean simplifies to \[\begin{align*} \lim_{T \to \infty} T\mathsf{Va}\left(\overline{Y}_T\right) = \tau^2 \left\{1+2\sum_{t=1}^\infty \rho_t\right\}. \end{align*}\] which is a function of

  • the unconditional variance \(\tau^2\)
  • the lag-\(h\) autocorrelation, \(\mathsf{Cor}(Y_{t}, Y_{t+h})=\rho_h.\)

Correlogram

For a weakly stationary process, it is instructive to look at the autocorrelation \(\rho_h=\mathsf{Cor}(Y_t, Y_{t+h}).\)

Figure 3: Correlogram of two Markov chains. These plots, often called acf or autocorrelation functions, show the lag-h sample autocorrelation against lag number.

Effective sample size

The effective sample size is, loosely speaking, the equivalent number of observations if the \(B\) marginal posterior draws where independent and more formally \[ \mathsf{ESS} = \frac{B}{\sum_{t=-\infty}^\infty \rho_t}= \frac{B}{\left\{1+2\sum_{t=1}^\infty \rho_t\right\}} \] where \(\rho_t\) is the lag \(t\) correlation, since \(\rho_{t} = \rho_{-t}\) and \(\rho_0=1.\)

Effective sample size for AR(1)

The lag-\(h\) correlation of the stationary autoregressive process of order 1 is \(\phi^h\), so \[1+2\sum_{h=1}^\infty \rho_h = 1 + 2 \left(\frac{1}{1-\phi}-1\right) = \frac{1+\phi}{1-\phi}.\]

whereas the unconditional variance of the \(\mathsf{AR}(1)\) process is

\[\tau^2 = \frac{\sigma^2}{1-\phi^2}.\]

Inefficiency curve for AR(1)

Figure 4: Left: scaled asymptotic variance of the sample mean for AR(1) (full line) and independent observations with unit marginal variance (dashed). Right: variance ratio for positive correlations for selected range.

To get the same precision for the mean of \(\mathsf{AR}(1)\) process with \(\phi \approx 0.75\) than with i.i.d. data, we would need 7 times as many observations.

Morale of the story

The price to pay for having correlated samples is

inefficiency

The higher the autocorrelation, the larger the variability of our estimators.

Standard error of mean for first-order autoregressive process

phi <- 0.8 # lag-one autocorrelation
sigma <- 2 # std. deviation of white noise
B <- 1e4L # chain length
x <- arima.sim( # generate AR(1)
  model = list(ar = phi),
  n = B,
  rand.gen = rnorm,
  sd = sigma
)
# Effective sample size / B
ess_factor_ar1 <- 1 / (-1 + 2 / (1 - phi))
# Unconditional variance tau_sq, obtained from var(x)
uvar_ar1 <- sigma^2 / (1 - phi^2)
# Theoretical std. error of sample mean
se_mean_ar1 <- sqrt(uvar_ar1 / B / ess_factor_ar1)

Standard errors estimation via nonparametric spectral analysis

Suppose that \(X_1, \ldots, X_B\) is a mean-stationary sequence.

Heidelberger & Welch (1981) recall that the variance of the sample mean is \(S(0) / B,\) where \(S(\cdot)\) is the spectrum evaluated at zero and \[ S(0) = \sum_{h=-\infty}^\infty \gamma_h. \] They propose to approximate \(S(0)\) by

  • calculating the raw periodogram (unsmoothed)
  • fitting a gamma GLM to the log periodogram with a polynomial trend.

Standard errors estimation via parametric spectral analysis

Any continuous spectral density function can be approximated by an \(\mathsf{AR}(p)\) process (for \(p\) potentially infinitely large).

For an autoregressive process of order \(p\), the spectral density of the \(\mathsf{AR}(p)\) is (Percival & Walden, 2020, eq. 446b) \[ S(0) = \frac{\sigma^2}{\left(1 - \sum_{j=1}^p \phi_j\right)^2}. \] This suggests

  • Fitting an \(\mathsf{AR}(p)\) model (selecting \(p\) via \(\mathsf{AIC}\)) using Yule–Walker or any other method.
  • Calculating the spectral density at zero with the plug-ins \(\widehat{\phi}_1, \ldots, \widehat{\phi}_p\) and \(\widehat{\sigma}^2.\)

Spectral estimation in R

# Heidelberg and Welch (1981)
coda::spectrum0(x)
$spec
[1] 83.04257
# Code from coda::effectiveSize
ar_mod <- ar.yw(x)
suvar <- ar_mod$var.pred / (1 - sum(ar_mod$ar))^2
se3 <- sqrt(suvar / B) # variance of sample mean of AR(p)
ess3 <- var(x) / suvar * B

Standard errors for posterior means using batch means

Geyer (2011) recommends to segment the time series into batches

  1. Break the chain of length \(B\) (after burn in) in \(K\) blocks of size \(\approx K/B.\)
  2. Compute the sample mean of each segment.
  3. Compute the standard deviation of the segments mean.
  4. Rescale by \(K^{-1/2}\) to get standard error of the global mean.

Illustration of batch means

Figure 5: Calculation of the standard error of the posterior mean using the batch method.

Convergence of Markov chains

If a Markov chain is irreducible and acyclic, it has a unique stationary distribution.

  • irreducibility: means that the chain can move from anywhere to anywhere, so it doesn’t get stuck in part of the space forever.
  • acyclic: cyclical chains loop around and visit periodically a state

The ergodic theorem holds even if the chain is cyclical.

Examples of cyclical or reducible chains

Consider discrete Markov chains over the integers \(1, 2, 3\) with transition matrices

\[ P_1 = \begin{pmatrix} 0.5 & 0.3 & 0.2 \\ 0 & 0.4 & 0.6 \\ 0 & 0.5 & 0.5 \end{pmatrix}, \quad P_2 = \begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix}. \] Chain 1 is reducible to \(\{2, 3\}\), chain 2 is cyclical.

Convergence and stationary distribution

Consider a transition \(P\) on \(1, \ldots, 5\) defined as \[ P = \begin{pmatrix} \frac{2}{3} & \frac{1}{3} & 0 & 0 & 0 \\ \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 & 0 \\ 0 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 \\ 0 & 0 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \\ 0 & 0 & 0 & \frac{1}{3} & \frac{2}{3} \\ \end{pmatrix} \] The stationary distribution is the value of the row vector \(\boldsymbol{p},\) such that \(\boldsymbol{p} = \boldsymbol{p}\mathbf{P}\) for transition matrix \(\mathbf{P}\): we get \(\boldsymbol{p}=(1,2,2,2,1)/8.\)

Convergence of Markov chains

Figure 6: Discrete Markov chain on integers from 1 to 5, with traceplot of 1000 first iterations (left) and running mean plots of sample proportion of each state visited (right).

Summary

  • Monte Carlo methods often rely on easy-to-simulate distributions \(q(x)\) to draw from a target \(p(x).\)
  • Rejection sampling: find

References

Botev, Z., & L’Écuyer, P. (2017). Simulation from the normal distribution truncated to an interval in the tail. Proceedings of the 10th EAI International Conference on Performance Evaluation Methodologies and Tools on 10th EAI International Conference on Performance Evaluation Methodologies and Tools, 23–29. https://doi.org/10.4108/eai.25-10-2016.2266879
Devroye, L. (1986). Non-uniform random variate generation. Springer. http://www.nrbook.com/devroye/
Geyer, C. J. (2011). Introduction to Markov chain Monte Carlo. In S. Brooks, A. Gelman, G. Jones, & X. L. Meng (Eds.), Handbook of Markov chain Monte Carlo (pp. 3–48). CRC Press. https://doi.org/10.1201/b10905
Heidelberger, P., & Welch, P. D. (1981). A spectral method for confidence interval generation and run length control in simulations. Communications of the ACM, 24(4), 233–245. https://doi.org/10.1145/358598.358630
Percival, D. B., & Walden, A. T. (2020). Spectral analysis for univariate time series. Cambridge University Press. https://doi.org/10.1017/9781139235723