1 Introduction and background

1.1 Motivation

The first part of the module is focused on random sampling and Monte Carlo (MC) methods. Loosely speaking, the MC methods use random sampling to make numerical estimates. It is widely used in various fields to simulate complex systems and predict outcomes when exact solutions are difficult or impossible to find. By running a large number of random trials, the method helps in making informed decisions based on probabilities. It is particularly valuable in areas like finance, physics, and engineering.

While MC methods can be used to solve a wide range of problems. Its key motivation can be summarised as the following expectation evaluation problem. Let $x$ be a random variable defined on $\mathbb R^d$ , following distribution $p(x)$ , and let $f(\cdot)$ be a real-valued function defined on $\mathbb R^d$ . We want to compute the expectation of $f(x)$ (clearly $f(x)$ is a random variable defined on $\mathbb R$ ): $\mathcal{I}=\mathbb E[f]=\int f(x)p(x)dx.\tag{1.1}$ Mathematically it can be understood as an integral evaluation, and we should note that, in all but the simplest situations, this integral does not admits an analytical form and needs to be evaluated numerically. Conventional numerical integration method can not efficiently evaluate the integral when the dimensionality $d$ is high. MC methods are especially suited for solving the expectation evaluation problem.

Original page 02

1.2 Review of some Statistics Basics

We first review some basics in Statistics that will be used in this module.

Random sample: given a distribution $F$ , let $X_1,\ldots,X_n$ be $n$ independent random variables following the distribution $F$ , and we call $X_1,\ldots,X_n$ a (independent and identically distributed, i.i.d.) random sample of size $n$ from distribution $F$ . $n$ is often called the sample size.

Loosely speaking, $X_1,\ldots,X_n$ can be understood as a set of realizations or data points following the distribution $F$ .

Parameter: a parameter of a distribution $F$ is a value that describes certain characteristic of the distribution $F$ , such as the distribution mean or variance.

Statistic: a statistic is a function of the sample $X_1,\ldots,X_n$ . E.g. the sample mean is a statistic: $\hat X=\frac1n\sum_{i=1}^nX_i.$ Estimator: let $X_1,\ldots,X_n$ be a sample from distribution $F$ , and $\theta$ is a parameter of $F$ . A (point) estimator of $\theta$ is a function (or a statistic) of the sample $X_1,\ldots,X_n$ : $\hat\theta=g(X_1,\ldots X_n),$ that is usually used to estimate parameter $\theta$ .

When we have an estimator, it is important to know how accurate it is, which is usually measure by the so-called mean-square error (MSE). Suppose that we have a sample $X_1,\ldots,X_n$ and an estimator $\hat\theta$ of parameter $\theta$ . First it should be clear that since $\hat\theta$ is a function of $X_1,\ldots,X_n$ , it is a random variable. The MSE is defined as, $\mathop{\mathrm{MSE}}[\hat\theta]=\mathbb E_{\hat\theta}[(\hat\theta-\theta)^2].\tag{1.2}$ From Eq. (1.2) it is easy to derive that $\begin{aligned} \mathop{\mathrm{MSE}}[\hat\theta]&=\mathbb E[(\hat\theta-\theta)^2]=\mathbb E[(\hat\theta-\mathbb E\hat\theta+\mathbb E\hat\theta-\theta)^2]\\ &=\mathbb E[(\hat\theta-\mathbb E\hat\theta)^2]+(\mathbb E\hat\theta-\theta)^2\\ &=\mathop{\mathrm{Var}}[\hat\theta]+(\mathbb E\hat\theta-\theta)^2 \end{aligned}\tag{1.3}$

Original page 03

It can be seen here that the MSE consists of two parts: variance and bias. Bias of the estimator $\hat\theta$ is the difference between this estimator’s expected value and the true value of the parameter being estimated: $\mathop{\mathrm{Bias}}[\hat\theta]=\mathbb E\hat\theta-\theta$ .

Variance of the estimator $\mathop{\mathrm{Var}}[\hat\theta]=\mathbb E[(\hat\theta-\mathbb E\hat\theta)^2$ , loosely speaking measures how much uncertainty in the estimator.

Unbiased estimator: an estimator $\hat\theta$ of parameter $\theta$ is said to be unbiased if $\mathop{\mathrm{Bias}}[\hat\theta]=0$ .

2 Basic Monte Carlo Method

Problem statement: let $x$ be a random variable defined on $\mathbb R^d$ , following distribution $p(x)$ , and let $f(\cdot)$ be a real-valued function defined on $\mathbb R^d$ . We want to compute the expectation of $y=f(x)$ (clearly $y$ is also a random variable): $\mathcal{I}=\mathbb E[y]=\int f(x)p(x)dx.\tag{2.1}$ How to compute $\mathcal{I}$ ? The problem can be understood as an integration problem and so numerical integration techniques can be used in principle. But here we want to use statistical methods.

The idea is the following: suppose we have a sample of $y$ : $y_1,\ldots,y_n$ , and we can estimate $\mathcal{I}$ as $\mathcal{I}\approx\widehat{\mathcal{I}}=\frac1n\sum_{i=1}^ny_i.\tag{2.2}$ However, in practice, we do not have such a sample, what should we do? We generate the sample using a computer!

The Monte Carlo (MC) method: first generate a $n$ -size sample of random variable $x$ : $x_1,\ldots,x_n$ , for each $x_i$ we can compute $f(x_i)$ , and estimate $\mathcal{I}$ as $\mathcal{I}\approx\widehat{\mathcal{I}}=\frac1n\sum_{i=1}^nf(x_i).\tag{2.3}$ This method is called the MC method and $\widehat{\mathcal{I}}$ is called the MC estimator of the expectation $\mathcal{I}$ . In what follows we will for simplicity omit the output

Original page 04

variable $y$ and only write the MC method in terms of the function $f$ : $\mathcal{I}\approx\widehat{\mathcal{I}}=\frac1n\sum_{i=1}^nf(x_i).\tag{2.4}$ Next we shall examine how accurate the MC estimator is, which is measured by the MSE (recall that it is made up of bias and variance). First we will show that the MC estimator is unbiased: $\begin{aligned} \mathbb E(\widehat{\mathcal{I}})&=\mathbb E\left(\frac1n\sum_{i=1}^nf(x_i)\right)=\frac1n\sum_{i=1}^n\mathbb E(f(x_i))\\ &=\mathbb E(f)=\mathcal{I}. \end{aligned}\tag{2.5}$ This means that the MSE of the MC estimator only depends on its variance. Next we look at the variance of MC estimator: $\begin{aligned}\mathop{\mathrm{Var}}[\widehat{\mathcal{I}}]&=\mathop{\mathrm{Var}}\left[\frac1n\sum_{i=1}^nf(x_i)\right]=\frac1{n^2}\sum_{i=1}^n\mathop{\mathrm{Var}}[f(x_i)]\\ &=\frac1{n^2}n\mathop{\mathrm{Var}}[f]=\frac{\mathop{\mathrm{Var}}[f]}n, \end{aligned}$ and it follows that the standard deviation of $\widehat{\mathcal{I}}$ is $\frac{\mathop{\mathrm{Std}}[f]}{\sqrt n}$ , which means that the standard deviation decreases at the rate of $\sqrt n$ . Thus people often say that the MC method has a convergence rate of $1/\sqrt n$ .

Now we consider the question: how to improve the accuracy of a MC estimator? Since MC estimator is unbiased, we only need to reduce its variance. Since $\mathop{\mathrm{Var}}[f]$ is fixed, the only way to reduce the variance is to increase $n$ , the sample size. For example, if we want to reduce the standard deviation to its $1/10$ , we need to use 100 times of the samples. Recall that each sample is actually computed via computer, and so that larger samples size simply means higher computational cost. Now an important question is can we reduce the variance without increasing the sample size? This can be done with the so-called variance reduction techniques.

Exercise: use MC method to evaluate the integral $\int_0^1x^2dx$ , $\int_0^1x^3dx$ and $\int_0^1\sin(\pi x)dx$ . Try different sample sizes $n=10,100,1000$ and record the estimation error (the difference between your estimated value and the truth). Hint: the first integral can be written as $\int x^2p(x)dx$ where $p(x)$ is a uniform distribution $U[0,1]$ and the other two are similar.

Original page 05

3 Importance Sampling

3.1 Importance sampling

Recall that the variance of a MC estimator is $\mathop{\mathrm{Var}}[\widehat{\mathcal{I}}]=\frac{\mathop{\mathrm{Var}}[f]}n,$ and if we want to reduce it we can increase the denominator or the sample size $n$ . However, a larger sample size means higher computational cost, and to this end it is natural to ask can we reduce the variance by decreasing the numerator $\mathop{\mathrm{Var}}[\widehat{\mathcal{I}}]$ ? But remember that we have said for a given problem $\mathop{\mathrm{Var}}[\widehat{\mathcal{I}}]$ is fixed, and so how can we change it? This can be done by using the variance reduction techniques, which gives us the flexibility to change $\mathop{\mathrm{Var}}[\widehat{\mathcal{I}}]$ . In this section we shall introduce one of the most popular variance reduction technique – importance sampling (IS).

The idea of IS is quite simple. Recall that the goal of MC is to compute the expectation $\mathbb E[f]=\int f(x)p(x)dx$ , which can be rewritten as $\begin{aligned} \mathbb E[f]&=\int f(x)p(x)dx=\int f(x)\frac{p(x)}{q(x)}q(x)dx\\ &=\int f(x)w(x)q(x)dx=\mathbb E_q[fw], \end{aligned}$ where $w(x)=\frac{p(x)}{q(x)}$ and $q(x)$ is an arbitrary distribution of $x$ . It should be clear that $\mathbb E(f)=\mathbb E_q(fw)$ and next we shall apply MC to $\mathbb E_q(fw)$ , yielding $\widehat{\mathcal{I}}_{IS}=\frac1n\sum_{i=1}^nf(x_i)w(x_i),$ where the samples $x_1,\ldots,x_n$ are drawn according to the new distribution $q(x)$ . Here $\widehat{\mathcal{I}}_{IS}$ is called the IS estimator of $\mathbb E[f]$ , and $q(x)$ is called the IS distribution, the biasing distribution, the alternative distribution or the proposal distribution in literature. It is trivial to show that it is an unbiased estimator for $\mathbb E[f]$ : $\mathbb E[\widehat{\mathcal{I}}_{IS}]=\mathbb E_q[fw]=\mathbb E[f].$

Original page 06

Next we look at the variance of IS estimator: $\begin{aligned} \mathop{\mathrm{Var}}[\widehat{\mathcal{I}}_{IS}]&=\mathop{\mathrm{Var}}\left[\frac1n\sum_{i=1}^nf(x_i)w(x_i)\right]=\frac1n\mathop{\mathrm{Var}}[fw]\\ &=\frac1n\bigl(\mathbb E_q[f^2w^2]-(\mathbb E_q[fw])^2\bigr)\\ &=\frac1n\left(\int f^2(x)w^2(x)q(x)dx-\mathbb E^2[f]\right)\\ &=\frac1n\left(\int f^2(x)w(x)p(x)dx-\mathbb E^2[f]\right) =\frac1n(\mathbb E[f^2w]-\mathbb E^2[f]), \end{aligned}\tag{3.1}$ which means that we now have the flexibility to reduce the variance by choosing an appropriate distribution $q(x)$ . To show this we can consider the best case scenario: there exists an optimal IS distribution $q(x)$ which may result in zero variance (that is, no error at all!). Let us verify that if $q(x)=\frac{f(x)p(x)}{\int f(x)p(x)dx},$ we have $w(x)=\frac{p(x)}{q(x)}=\frac{\int f(x)p(x)dx}{f(x)},$ and the variance of the resulting IS estimator is, $\mathop{\mathrm{Var}}[\widehat{\mathcal{I}}_{IS}]=\frac1n(\mathbb E[f^2w]-\mathbb E^2[f])=\frac1n(\mathbb E^2[f]-\mathbb E^2[f])=0,$ which means that as long as $n>0$ , the estimator is exact! But we note here that this optimal IS distribution can not be used directly, as it actually requires the knowledge of the sought after expectation $\mathbb E[f]$ . Nevertheless it can still be used to help design efficient IS estimators (to be discussed later).

In summary we can see that choosing a good IS distribution is the key of the IS method. However making a good choice of the IS distribution can be extremely hard: sometimes it is problem-dependent and so we have to design an IS distribution by analyzing the underlying problem, and sometimes it can be obtained using some rather generic schemes.

3.2 Small Probability Estimation

In this section we consider a very special but important type of problems where MC and IS methods can apply – probability estimation. That is, we

Original page 07

want to estimate the probability of a random event $P(x\in A)$ , and we can rewrite it in the form of an expectation by defining an indicator function: $P(x\in A)=\int_{x\in A}p(x)dx=\int I_A(x)p(x)dx=\mathbb E[I_A(x)]\tag{3.2}$ where $I_A(x)=\begin{cases}1&x\in A,\\0&x\notin A.\end{cases}$ Now we can use the MC method to compute the probability of interest: $\hat P=\frac1n\sum_{i=1}^n I_A(x_i),$ where $x_1,\ldots,x_n$ are drawn from the distribution $p(x)$ .

We shall provide two examples here.

Example: coin flips.

Experiment: 10 coin flips. Question: What is the probability of 9 heads?

Answer: it follows a binomial distribution and can be calculated directly: $\approx0.01$

But, how do we compute it with experiments?

Solution: Use an unfair coin!

Optimal: Use a coin that gives heads 90% of the time.

Correct for biased coin by using likelihood ratio: If on a flip one gets heads, multiply by $0.5/0.9$

If on a flip one gets tails, multiply by $0.5/0.1$

This corrects statistics: get results for a fair coin

In terms of MC/IS, the original distribution is Bernoulli with the probability of 1 (corresponding to head) being $p=0.5$ , and the IS distribution is also Bernoulli with $p=0.9$ .

Example: $p(x_1,x_2)=N(0,I_2)$ and $A=\{x|(x_1-1)^2+(x_2-1)^2<1\}$ . Calculate $P(x\in A)$ . Use MC to calculate the probability in the examples above.

What if the region $A$ becomes $A=\{x|(x_1-5)^2+(x_2-5)^2<1\}$ ? Why it is rare? (Use a graph to show that majority of the samples do not fall in the region $A$ .) Use IS (try different parameter values) to calculate the

Original page 08

probability. Hint: let the IS distribution to be Gaussian whose mean is to be determined. This technique is called mean-shift, and the key question is where we want to shift the distribution? Answer: we want to encourage the samples to fall in the important region. (Use graph to explain the "important region".)

3.3 Variance minimisation and exponential tilting

3.3.1 Variance minimisation

As has been mentioned earlier, it is highly useful to have some generic schemes for finding a good IS distribution. Here we will introduce a very straightforward idea – to find an IS distribution that leads to small estimation variance, a strategy often referred to as variance minimisation. Rather than consider all possible choices for the sampling distribution $q(x)$ , one strategy is to restrict the set of $q(x)$ we consider to some parametrised family of distributions and minimize the estimator variance over this family. So we assume we have a family of distributions $q(x,\theta)$ where $\theta$ parametrises the family. Both $x$ and $\theta$ are multi-dimensional. The weighting function is $w(x,\theta)=\frac{p(x)}{q(x,\theta)}$ Recall that the variance of the IS estimator is, $\sigma^2(\theta)=\int f(x)^2w(x,\theta)p(x)dx-\mathbb E[f]^2$ As mentioned earlier, we can find a good IS distribution by minimising this as a function of $\theta$ (note that it is equivalent to minimizing the first integral on the right-hand-side of the equation), namely, $\min_{\theta\in\Theta}F(\theta):=\int f(x)^2w(x,\theta)p(x)dx\tag{3.3}$ where $\Omega$ is a prescribed feasible set of $\theta$ . Solving the optimisation problem directly is usually difficult, largely due to the integral in the objective function. In practice, we often solve the MC estimation of the original objective function: $\min_{\theta\in\Theta}\hat F(\theta):=\frac1n\sum_{i=1}^nf(x_i)^2w(x_i,\theta)\tag{3.4}$

Original page 09

where $x_1,\ldots,x_n$ are drawn from $p(x)$ . Next we need to discuss how to design the parametrised distribution family.

3.3.2 Exponential tilting

Exponential tilting is a technique used construct importance sampling distribution. Specifically the method uses the so-called tilted distribution $q(x;\theta)$ defined as: $q(x;\theta)=\frac{p(x)e^{\theta^Tx}}{c(\theta)}\tag{3.5}$ where $\theta$ is called the tilting parameter, and $c(\theta)=\mathbb E[e^{\theta^TX}]$ is the normalizing constant that ensures $q(x;\theta)$ is a valid probability distribution. It should be clear that in this case, $\theta$ and $x$ are of the dimensionality.

By sampling from the tilted distribution $q(x;\theta)$ , the importance sampling estimator is: $\begin{aligned} \mathbb E[f(X)]&=\mathbb E_q\left[f(X)\frac{p(X)}{q(X;\theta)}\right] =\mathbb E_q\left[f(X)\frac{c(\theta)]}{e^{\theta^TX}}\right]\\ &\approx\frac1n\sum_{i=1}^nf(x_i)w(x_i,\theta) \end{aligned}$ where $w(x)=c(\theta)/e^{\theta^Tx}$ . The optimal value of $\theta$ can be determined by using the variance minimisation method, i.e., by substituting Eq. (3.5) into Eq. (3.3) and solving for $\theta$ . A difficulty here is that, for an arbitrary distribution $p(x)$ , the tilted distribution $q(x)$ may not be a standard (commonly used) distribution, and thus sampling from it can be challenging.

3.3.3 Exponential family

Definition 3.1. A parametrised distribution family in the form of $p_{Exp}(x,\theta)=h(x)\exp(\eta(\theta)^TT(x)-A(\theta)),\tag{3.6}$ where $h(x)$ and $A(\theta)$ are scalar-valued functions, and $\eta(\theta)$ and $T(x)$ are vector-valued functions, is called an exponential family. As a special case of the exponential family, when $\eta(\theta)=\theta$ and $T(x)=x$ , it is called a natural exponential family.

Theorem 1. A tilted distribution is in a natural exponential family.

Original page 10

Proof: A tilted distribution can be written as, $q(x;\theta)=\frac{p(x)e^{\theta^Tx}}{c(\theta)}=p(x)\exp[\theta^Tx-\log c(\theta)]\tag{3.7}$ and that is $q(x,\theta)$ is a natural exponential family with $h(x)=p(x)$ and $A(\theta)=\log c(\theta)$ .

Theorem 2. If the original distribution $p(x)$ is in a natural exponential family, its tiled distribution is in the same natural exponential famility.

Proof: We have $p(x)=h(x)\exp(\theta_0^Tx-A(\theta_0))$ and it should be clear that $A(\theta)=\log\int h(x)\exp(\theta^Tx)dx$ Its tilted distribution is now $\begin{aligned} q(x;\theta)&=\frac{p(x)e^{\theta^Tx}}{c(\theta)} =h(x)\exp[\theta_0^Tx-A(\theta_0)]\frac{e^{\theta^Tx}}{c(\theta)}\\ &=h(x)\exp[(\theta_0+\theta)^Tx-B(\theta,\theta_0),] \end{aligned}\tag{3.8}$ where $B(\theta,\theta_0)$ is the normalisation constant. Next we determine the value of $B$ : $B(\theta,\theta_0)=\int h(x)\exp[(\theta_0+\theta)^T]xdx=A(\theta_0+\theta).$ Thus we show that $q(x,\theta)$ is from the same natural exponential family of $p(x)$ .

Finally a few remarks:

If the original distribution is in a common natural exponential family, the exponential tilted method is from the same family and so is easy to use.
While it is a special case, the natural exponential families are indeed widely used in practice, which makes the method (Variance Minimisation+Exponential Tilting+Natural Exponential Family) highly useful.
That said, the variance minimisation method is often computationally intensive. In Section 4, we will present an alternative method (the Cross Entropy method), which is often more efficient to implement.

Original page 11

3.4 Self-normalised importance sampling

Self-normalized importance sampling (SNIS) is a method used when the target density $p(x)$ and/or the proposal distribution $q(x)$ are only known up to a constant factor. This is common in many practical problems, such as Bayesian inference. The idea behind SNIS is to avoid needing to compute the normalizing constant for $p(x)$ , which can be unknown or computationally intractable.

Recall that our goal is to estimate the integral $\mathbb E[f]=\int f(x)p(x)dx.$ When $p(x)$ is known up to a constant, i.e., $p(x)=c_pp_0(x)$ , and $q(x)=c_qq_0(x)$ is the proposal distribution, the self-normalized importance sampling estimator is constructed as follows:

Rewrite the expectation: $\begin{aligned} \mathbb E_p[f(x)]&=\frac{\int f(x)\frac{p(x)}{q(x)}q(x)dx}{\int\frac{p(x)}{q(x)}q(x)dx} =\frac{\int f(x)\frac{p_0(x)}{q_0(x)}q(x)dx}{\int\frac{p_0(x)}{q_0(x)}q(x)dx}\\ &=\frac{\mathbb E_q[f(x)w(x)]}{\mathbb E_q[w(x)]} \end{aligned}$ where $w(x)=\frac{p_0(x)}{q_0(x)}$ .
The SNIS estimator is given by: $\widehat{\mathcal{I}}_{SNIS}=\frac{\sum_{n=1}^Nf(x_n)w(x_n)}{\sum_{i=1}^Nw(x_n)}$ where $x_i$ are the samples drawn from the proposal distribution $q(x)$ .

Next we will look at the MSE of the SNIS estimator. Unlike the standard IS estimator, it is biased. This can be seen by considering the special case $n=1$ : $\mathbb E_q[\widehat{\mathcal{I}}_{SNIS}]=\mathbb E_q\left[\frac{f(x_1)w(x_1)}{w(x_1)}\right] =\int f(x)q(x)dx\ne\int f(x)p(x)dx.$

Original page 12

Definition 3.2. An estimator $\hat\theta_n$ of a parameter $\theta$ is said to be asymptotically unbiased if its expected value converges to the true parameter as the sample size $n$ approaches infinity. That is, $\lim_{n\to\infty}\mathbb E[\hat\theta_n]=\theta.$ We can show that the estimator $\widehat{\mathcal{I}}_{SNIS}$ is an asymptotically unbiased estimator of $\mathbb E[f]$ . That is to say, although the SNIS estimator is biased, the bias decreases as the number of samples grows.

Proof. $\widehat{\mathcal{I}}_{SNIS}= \frac{\frac1N\sum_{n=1}^Nf(x_n)\frac{p(x_n)}{q(x_n)}}{\frac1N\sum_{n=1}^N\frac{p(x_n)}{q(x_n)}}$ Now apply the strong law of large numbers to both numerator and denominator separately. The numerator converges to $\mathbb E[f]$ and the denominator converges to 1, almost surely. It follows immediately that $\mathbb E[\widehat{\mathcal{I}}_{SNIS}]$ converges to $\mathbb E[f]$ as $N$ increases.

In case your need, Strong Law of Large Numbers (SLLN) says:

Theorem 3. Let $\{X_n\}$ be a sequence of independent and identically distributed (i.i.d.) random variables with common mean $\mu=\mathbb E[X_n]$ and finite variance. Then the sample average converges almost surely to the expected value $\mu$ , i.e., $\frac1n\sum_{i=1}^nX_i\xrightarrow{\mathrm{a.s.}}\mu\quad\text{as}\quad n\to\infty.$ This means that $P\left(\lim_{n\to\infty}\frac1n\sum_{i=1}^nX_i=\mu\right)=1.$

Original page 13

4 The cross entropy method

4.1 The motivating idea

As is mentioned above, it is also possible to derive a good IS distribution using some generic approaches. The cross entropy (CE) method is one of them. Recall that the the IS estimator admits an optimal IS distribution, which can yield zero variance: $q^*(x)=\frac{f(x)p(x)}{\mathbb E[f]},$ and we know that we can not directly use this optimal distribution as it depends on the sought-after $\mathbb E[f]$ . If we can find a distribution that is “close” to the optimal one, intuitively it should also be a rather good IS distribution. This is the basic idea behind the CE method. Now we need some distance measure that can quantify the similarity between two distributions.

4.2 The Kullback-Leibler divergence

A very popular distance for measuring the similarity between two distributions is the Kullback-Leibler divergence (KLD) for continuous random variables: $D(p_1,p_2)=\int\log\frac{p_1(x)}{p_2(x)}p_1(x)dx,\tag{4.1}$ which is also known as the cross entropy between the two distributions. The KLD for discrete random variables is $D(p_1,p_2)=\sum_{x_i\in S}\log\frac{p_1(x_i)}{p_2(x_i)}p_1(x_i)dx,\tag{4.2}$ where $p_1$ and $p_2$ are two probability mass functions defined on state space $S$ . In what follows we shall focus on the KLD for the continuous random variables, while noting that most of the conclusions apply to the discrete version as well.

Next we shall look at some important properties of KLD.

It is non-negative

Original page 14

it is zero only if $p_1=p_2$
however, it is not symmetric
it also does not satisfy the triangle property which requires $D(p_1,p_3)\leq D(p_1,p_2)+D(p_2,p_3)$ .

The last two properties show that KLD is not a well-defined "distance metric"; nevertheless KLD is widely used in measuring the distance between two distributions. Now we prove these properties of KLD.

Proof : (1) and (2): we have, $\begin{aligned} D(p_1,p_2)&=\int\log\frac{p_1(x)}{p_2(x)}p_1(x)dx\\ &=-\int\log\frac{p_2(x)}{p_1(x)}p_1(x)dx \geq-\log\int\frac{p_2(x)}{p_1(x)}p_1(x)dx=0, \end{aligned}\tag{4.3}$ which is a consequence of Jensen’s inequality. (Jensens’s inequality: if $g(x)$ is a convex function, and $x$ is a random variable, $\mathbb E[g(x)]\leq g(\mathbb E[x]])$ , and if $g(x)$ is strictly convex, the equality holds if and only if $x$ is a constant.) Because $\log x$ is a strictly convex function, the equality holds if and only if $p_1(x)/p_2(x)$ is a constant (this implies $p_1=p_2$ , why?).

(3) it is obvious.

(4) $D(p_1,p_3)$ can be $<$ , $=$ or $>$ than $D(p_1,p_2)+D(p_2,p_3)$ , which can be proved by example. We shall use the example of Bernoulli distribution. First let let $p_1$ be a Bernoulli distribution with parameter $a_1$ and $p_2$ be that with parameter $a_2$ . First we shall derive that the KLD between two Bernoulli distributions is $\begin{aligned} D(p_1,p_2)&=\sum_{x_i\in\{0,1\}}\log(p_1(x_i)/p_2(x_i))p_1(x_i)\\ &=\log(a_1/a_2)\times a_1+\log((1-a_1)/(1-a_2))\times(1-a_1). \end{aligned}\tag{4.4}$ Now consider three Bernoulli distributions: $p_1$ with $a_1=0.1$ , $p_2$ with $a_2=0.5$ and $p_3$ with $a_3=0.9$ . Using Eq. (4.4) we can compute $D(p_1,p_2)=0.3681$ , $D(p_2,p_3)=0.5108$ and $D(p_1,p_3)=1.7578$ , which clearly violates the triangle inequality.

Original page 15

4.3 The CE algorithm

Next we shall apply the KLD to identify a good IS distribution and that is we seek a distribution that is “close” to the optimal distribution with respect to KLD. More precisely we shall find a good IS distribution (within a distribution family $\mathcal Q$ ) by minimizing its distance to the optimal IS distribution: $\min_{q\in\mathcal Q}D(q^*,q)=\int\log\frac{q^*(x)}{q(x)}q^*(x)dx\tag{4.5}$ $\begin{aligned} D(q^*,q)&=\int\log\frac{q^*(x)}{q(x)}q^*(x)dx\\ &=\int\log q^*(x)q^*(x)dx-\int\log q(x)q^*(x)dx \end{aligned}\tag{4.6}$ note here that the first time is a constant, and so $\min D(q^*,q)$ is equivalent to $\max_{q\in\mathcal Q}\int\log q(x)q^*(x)dx=\int\log q(x)\frac{f(x)p(x)}{\mathbb E[f]}dx$ or simply $\max_{q\in\Omega}\int\log q(x)f(x)p(x)dx,\tag{4.7}$ where amazingly nothing is unknown. In principle we now can find an IS distribution by solving the optimisation problem (4.7). However, directly optimizing over a density function can be rather difficult and in practice one often choose $q$ from a parameterized family. That is we take $q$ to be in the form of $q(x)=q(x|\theta)$ where $\theta$ represents some distribution parameter, and then turn Eq. (4.7) into $\max_{\theta\in\Omega}F(\theta\in\Omega)=\int\log q(x,\theta)f(x)p(x)dx,\tag{4.8}$ which is much easier to solve. Here $\Omega$ is the state space of the parameter $\theta$ . Moreover, in reality, we can not optimize the expectation and we optimize its MC estimator instead: $\max_{\theta\in\Omega}\hat F(\theta)=\frac1n\sum_{i=1}^n\log q(x_i,\theta)f(x_i),\tag{4.9}$ where $x_1,\ldots,x_n$ are drawn from distribution $p(x)$ . To summarize, the CE method proceeds as follows:

draw $n$ samples $x_1,\ldots,x_n$ from distribution $p(x)$ ;

Original page 16

compute $\theta$ by solving the optimisation problem: $\max_{\theta\in\Omega}\hat F(\theta)=\frac1n\sum_{i=1}^n\log q(x_i,\theta)f(x_i),$
draw samples $x_1,\ldots,x_n$ from $q(x|\theta)$ ;
estimate $\widehat{\mathcal{I}}_{IS}=\sum_{i=1}^nf(x_i)w(x_i)$ where $w(x)=p(x)/q(x|\theta)$ .

In practices, however, the algorithm above may not yield satisfactory IS distribution due to inaccurate estimation in Eq. (4.9) (as it is only an estimator of Eq. (4.8) which we really want to optimize), and so we will present an adaptive version of the algorithm which can gradually improve the IS distribution obtained. Roughly speaking, the idea is that, once we have computed an IS distribution denoted as $q(x|\theta^*)$ , we use it to re-estimate Eq. (4.8) with an IS estimator: $\begin{gathered} \max_{\theta\in\Omega}\hat F(\theta)=\frac1n\sum_{i=1}^n\log q(x_i|\theta)f(x_i)w(x_i),\\ w(x)=\frac{p(x)}{q(x|\theta^*)}, \end{gathered}\tag{4.10}$ where samples $x_1,\ldots,x_n$ are drawn from $q(x|\theta^*)$ ; we then solve the optimisation problem (4.10) and get a new and ideally better IS distribution; we then repeat this procedure until some stoping criteria are satisfied. The algorithm proceeds as follows:

initialization: choose $\theta_0$ , and let $k=0$ ;
while (stopping conditions not satisfied) do
1. draw $n$ samples $x_1,\ldots,x_n$ from distribution $q(x|\theta_k)$ ;
2. compute $\theta_{k+1}$ by solving the optimisation problem $\begin{gathered} \max_{\theta\in\Omega}\hat F(\theta)=\frac1n\sum_{i=1}^n\log q(x_i|\theta)f(x_i)w(x_i),\\ w(x)=\frac{p(x)}{q(x|\theta_k)}, \end{gathered}\tag{4.11}$
3. $k=k+1$ ;
draw samples $x_1,\ldots,x_n$ from $q(x|\theta_{k-1})$ ;
estimate $\widehat{\mathcal{I}}_{IS}=\sum_{i=1}^nf(x_i)w(x_i)$ where $w(x)=p(x)/q(x|\theta_{k-1})$ .

Original page 17

Commonly used stoping criteria include: 1) the difference between $\theta_{k-1}$ and $\theta_k$ is smaller than a threshold; 2) stop when a maximum number of iterations is reached.

4.4 CE for Gaussian IS distributions

In general Eq. (4.9) needs to be solved numerically, but for some special choices of the IS distribution, this problem can be solved analytically. Let us take $q$ to be an independent multivariate Gaussian distribution: $q(x)=\prod_{j=1}^dp_j(x_j),\qquad p_j(x_j)=N(\mu_j,\sigma_j^2),\tag{4.12}$ where the mean $\mu$ and variance $\sigma^2$ are the parameters to be determined. Substituting this Gaussian distribution into Eq. (4.8) yields $\begin{aligned} F(\theta)&=\int\log q(x,\theta)f(x)p(x)dx,\\ &=-\int\sum_{j=1}^d\left(\frac{(x_j-\mu_j)^2}{2\sigma_j^2}+\log\sigma_j\right)f(x)p(x)dx\\ &=-\sum_{j=1}^d\int\left(\frac{(x_j-\mu_j)^2}{2\sigma_j^2}+\log\sigma_j\right)f(x)p(x)dx \end{aligned}\tag{4.13}$ The optimal solution can be found by setting the gradient of $F(\cdot)$ in Eq. (4.9) to zero; first we consider $\mu_i$ (for $i=1\ldots d$ ): $\frac{\partial F}{\partial\mu_i}=\int(x_i-\mu_i)f(x)p(x)dx=0,\tag{4.14}$ which implies that $\mu_i=\frac{\int x_if(x)p(x)dx}{\int f(x)p(x)dx}.\tag{4.15}$ Next we take the gradient of $\sigma_i$ (for $i=1\ldots d$ ): $\frac{\partial F}{\partial\sigma_i}=-\int\left(\frac{(x_i-\mu_i)^2}{2}(-2)\sigma_i^{-3}+\frac1\sigma\right)f(x)p(x)dx=0,\tag{4.16}$ and by some basic calculus we get $\sigma_i^2=\frac{\int(x_i-\mu_i)^2f(x)p(x)dx}{\int f(x)p(x)dx}.\tag{4.17}$

Original page 18

In practice both $\mu_i$ and $\sigma_i$ are estimate via MC simulation: $\mu_i=\frac{\sum_{j=1}^nx_{i,j}f(x_j)}{\sum_{j=1}^nf(x_j)},\tag{4.18}$ $\sigma_i^2=\frac{\sum_{j=1}^n(x_{i,j}-\mu_i)^2f(x_j)}{\sum_{j=1}^nf(x_j)},\tag{4.19}$ where $x_{i,j}$ indicates the $i$ -th component of the $j$ -th sample, and $x_1,\ldots,x_n$ are samples drawn from $p(x)$ . As is discussed earlier, they can also be estimated with IS (let $q(x|\theta^*)$ be the present IS distribution): $\mu_j=\frac{\sum_{i=1}^nx_{i,j}f(x_{i,j})w(x_i)}{\sum_{i=1}^nf(x_i)w(x_i)},\tag{4.20}$ $\sigma_j^2=\frac{\sum_{i=1}^n(x_{i,j}-\mu_{i,j})^2f(x_i)w(x_i)}{\sum_{i=1}^nf(x_i)w(x_i)},\tag{4.21}$ where $w(x)=p(x)/q(x|\theta^*)$ and $x_1,\ldots,x_n$ are drawn from $q(x|\theta^*)$ .

4.5 Applications to rare probability estimation.

The algorithm above describes the CE method in the generic setting. Next we shall look at how it is applied to the problem of rare probability estimation. In principle the CE method can apply directly to this case, except that the function $f$ is an indicator function. However, an implementation issue is that, when the sought probability is extremely small, unless an unreasonably large number of samples are used, we may not see any samples fall in the region of interest (the support of the indicator function). To address the issue, we next present a multi-stage CE method. The idea is to construct a sequence of events with decaying probability which eventually arrive at the event of interest.

First for convenience’s sake, we formulate the event of interest as follows. Let $x$ be a random variable following distribution $p(x)$ , and $g(x)$ is a scalar function. We are interested in the probability $P(x|g(x)>\gamma)$ for a prescribed scalar $\gamma$ . As such the indicator function can be written as $I_\gamma(x)=\begin{cases}1,&g(x)\geq\gamma\\0,&g(x)<\gamma.\end{cases}\tag{4.22}$

Original page 19

Thus the probability of interest can be expressed as $P_\gamma=\mathbb E[I_\gamma]=\int I_\gamma(x)p(x)dx.$ As is mentioned earlier, if the probability of interest is too small, the standard CE will fail to compute an IS distribution. We use the following multi-stage CE method to address the issue. First of all, for a given sequence of $\gamma_1\leq\gamma_2\leq\ldots$ , we define the corresponding indicator function $I_{\gamma_k}(x)=\begin{cases}1,&g(x)\geq\gamma_k\\0,&g(x)<\gamma_k.\end{cases}\tag{4.23}$ Starting from an initial parameter $\theta_0$ and a not so small parameter $\rho$ (for example $\rho=0.05$ ), the multi-stage CE proceeds as follows,

initialization: choose $\theta_0,\rho$ ; let $k=1$ and $\gamma_0=+\infty$ ;
while $\gamma_{k-1}<\gamma$ do
- update $\gamma_k$ :
  1. generate $x_1,\ldots,x_n$ from $q(x|\theta_{k-1})$ ;
  2. compute $g(x_1),\ldots,g(x_n)$ ;
  3. let $\gamma_k$ the sample $(1-\rho)$ -th quantile of $\{x_1,\ldots,x_n\}$ ;
  4. if $\gamma_k>\gamma$ let $\gamma_k=\gamma$ ;
- update $\theta_k$ :
  - compute $\theta_{k+1}$ as the solution of the optimisation problem $\max_{\theta\in\Omega}\hat F(\theta)=\frac1n\sum_{i=1}^n\log q(x_i|\theta)I_{\gamma_k}(x_i)w(x_i),$ $w(x)=\frac{p(x)}{q(x|\theta_k)},$
- $k=k+1$ .

5 Other variance reduction methods

5.1 The control variates method

In this section, we discuss another popular variance reduction technique - the control variates method. Recall that our goal is to evaluate the expectation

Original page 20

$\mathbb E[f(x)]$ and now suppose we have another function $g(x)$ whose expectation is known in advance: $\mathbb E[g(x)]=g_0$ . We can now define a new function $f_{CV}(x)=f(x)+c(g(x)-g_0),$ where $c$ is an arbitrary constant. Next we can show that, $\mathbb E[f_{CV}(x)]=\mathbb E[f(x)+c(g(x)-g_0)]=\mathbb E[f(x)],$ and based on this we can construct the CV estimator of $\mathbb E[f]$ : $\widehat{\mathcal{I}}_{CV}=\frac1n\sum_{i=1}^n(f(x_i)+c(g(x_i)-g_0)),\tag{5.1}$ where $x_1,\ldots,x_n$ are drawn from $p(x)$ . It is easy to see that $\widehat{\mathcal{I}}_{CV}$ is an unbiased estimator of $\mathbb E[f]$ : $\mathbb E[\widehat{\mathcal{I}}_{CV}]=\mathbb E\left[\frac1n\sum_{i=1}^n(f(x_i)+c(g(x_i)-g_0))\right]=\mathbb E[f].$ The variance of $\widehat{\mathcal{I}}_{CV}$ is $\mathop{\mathrm{Var}}[\widehat{\mathcal{I}}_{CV}]=\mathop{\mathrm{Var}}\left[\frac1n\sum_{i=1}^nf_{CV}(x_i)\right]=\frac1n\mathop{\mathrm{Var}}[f_{CV}]$ where $\begin{aligned} \mathop{\mathrm{Var}}[f_{CV}]&=\mathop{\mathrm{Var}}[f+c(g-g_0)]=\mathop{\mathrm{Var}}[f+cg-cg_0]=\mathop{\mathrm{Var}}[f+cg]\\ &=\mathop{\mathrm{Var}}[f]+\mathop{\mathrm{Var}}[cg]+\mathop{\mathrm{Cov}}[f,cg]\\ &=\mathop{\mathrm{Var}}[f]+c^2\mathop{\mathrm{Var}}[g]+2c\mathop{\mathrm{Cov}}[f,g]. \end{aligned}\tag{5.2}$ It is easy to derive that the optimal value of $c$ is $c^*=-\frac{\mathop{\mathrm{Cov}}[f,g]}{\mathop{\mathrm{Var}}[g]},\tag{5.3}$ and the minimal variance is $\begin{aligned} \mathop{\mathrm{Var}}[f_{CV}]&=\mathop{\mathrm{Var}}[f]-\frac{\mathop{\mathrm{Cov}}^2[f,g]}{\mathop{\mathrm{Var}}[g]}\\ &=\mathop{\mathrm{Var}}[f]\left(1-\frac{\mathop{\mathrm{Cov}}^2[f,g]}{\mathop{\mathrm{Var}}[f]\mathop{\mathrm{Var}}[g]}\right) =\mathop{\mathrm{Var}}[f](1-\mathop{\mathrm{Corr}}^2[f,g]). \end{aligned}\tag{5.4}$

Original page 21

Now one can see that the stronger the correlation between $f$ and $g$ is, the smaller the estimator variance is. In other word, in the CV method, we shall seek a function $g$ that is as correlated to $f$ as possible, subject to the constraint that the expectation of $g$ is available.

Note that, in principle, the optimal value $c^*$ can be estimated using the Monte Carlo method: $\widehat{\mathcal{I}}_{MC}=\frac1n\sum_{i=1}^nf(x_i),\tag{5.5}$ $\hat a=\frac1n\sum_{i=1}^n(f(x_i)-\widehat{\mathcal{I}}_{MC})(g(x_i)-g_0),\tag{5.6}$ $\hat b=\frac1n\sum_{i=1}^n(g(x_i)-g_0)^2,\tag{5.7}$ $\hat c=\frac ab.\tag{5.8}$ Extension. The method can be extended to multiple control variate terms: namely suppose that we have a number of functions $g_l(x)$ with $\mathbb E[g_l]=g_0^l$ for $l=1,\ldots,L$ . We can define a function, $f_{CV}^L(x)=f(x)+\sum_{l=1}^Lc_l(g_l(x)-g_0^l),$ which then leads to a CV estimator: $\widehat{\mathcal{I}}_{CV}^L=\sum_{i=1}^n\left(f(x_i)+\sum_{l=1}^Lc_l(g_l(x_i)-g_0^l)\right)$ where $x_i$ are drawn from $p(x)$ . Next we $\begin{aligned} \mathop{\mathrm{Var}}[\widehat{\mathcal{I}}_{CV}^L]&=\mathop{\mathrm{Var}}\left[\frac1n\sum_{i=1}^nf_{CV}^L(x_i)\right]=\frac1n\mathop{\mathrm{Var}}[f_{CV}^L]\\ &=\frac1n\left(\mathop{\mathrm{Var}}[f]+\sum_{l=1}^L\sum_{l'=1}^L(c_lc_{l'})\mathop{\mathrm{Cov}}[g_l,g_{l'}] +2\sum_{l=1}^Lc_l\mathop{\mathrm{Cov}}[f,g_l]\right) \end{aligned}\tag{5.9}$ Taking the gradient of the equation above and setting it to zero, results in the following linear system (of $L$ equations), $\sum_{l'=1}^L\mathop{\mathrm{Cov}}[g_l,g_{l'}]c_{l'}=-\mathop{\mathrm{Cov}}[f,g_l],$

Original page 22

for $l=1,\ldots,L$ . Solving this linear system yields the optimal values of the coefficients $c_1,\ldots c_L$ .

Exercise: Use the control variates method to estimate the integral $\int_0^1\frac1{1+x}dx$ . (hint: treat the integral as an expectation with uniform distribution $U[0,1]$ ). Use $g(x)=1+x$ , $g_0=3/2$ (verify this yourself) and $c^*=0.4773$ . Compare your estimate against the true integral value (you can analytically evaluate it) and standard MC with the same amount of samples.

5.2 Antithetic variables

Definition 5.1 Random variables $X,Y$ on the same probability space are antithetic if they have the same distribution and their covariance is negative.

Suppose we want to compute $\mathbb E[f(x)]$ and we can find another function $g(x)$ such that $(f,g)$ is an antithetic pair. So $\mathbb E[g]$ is also equal to $\mathbb E[f]$ . Our Antithetic Monte Carlo (AMC) method is as follows. We generate $n$ independent samples $x_1,\ldots,x_n$ from $p(x)$ , and estimate $\mathbb E[f]$ as $\widehat{\mathcal{I}}_{AMC}=\frac1{2n}\sum_{i=1}^n(f(x_i)+g(x_i))\tag{5.10}$ Obviously, the mean of $\widehat{\mathcal{I}}_{AMC}$ is $\mathbb E[f]$ , so this is an unbiased estimator. Let $\rho$ denote the correlation of $f$ and $g$ . Since they have the same distribution, they have the same variance, denoted by $\sigma^2$ . The variance of our estimator is: $\begin{aligned} \operatorname{var}(\widehat{\mathcal{I}}_{AMC})&=\frac1{(2n)^2}n\operatorname{var}(f+g)\\ &=\frac1{2n}\sigma^2(1+\rho) \end{aligned}\tag{5.11}$ Obviously the variance of $\widehat{\mathcal{I}}_{AMC}$ is smaller than that of $\widehat{\mathcal{I}}_{MC}$ , as long as $\rho<0$ . Note however, there is additional cost of evaluating function $g(x)$ , and it is easy to see that if the cost of evaluating $g(x)$ is significantly higher that of evaluating $f(x)$ , the advantage of the method is limited. Therefore, the use of this method relies on the ability to find an antithetic function which is not computationally expensive.

Original page 23

5.3 Stratified sampling

Motivation

To motivate stratified sampling, consider the following simple example. We want to compute: $\mathcal{I}=\int_0^1f(x)\,dx\tag{5.12}$ On the interval $[0,1/2]$ , the function $f(x)$ is (nearly) constant, but on the interval $[1/2,1]$ the function varies significantly. Suppose we wanted to use Monte Carlo to compute the two integrals: $\int_0^{1/2}f(x)\,dx,\qquad\int_{1/2}^1f(x)\,dx\tag{5.13}$ The Monte Carlo for the first integral will have a much smaller variance ( $\approx0$ ) than the Monte Carlo for the second integral. Thus, it would make more sense to spend more time on the second integral, i.e., generate more samples for the second integral. However, under the usual Monte Carlo approach, we would sample from $U[0,1]$ and get approximately the same number of samples in $[0,1/2]$ and $[1/2,1]$ . The idea of stratified sampling is to divide the probability space into several regions and perform a Monte Carlo simulation for each region.

General setting

Let $\Omega$ be the probability space of distribution $p(x)$ . We partition it into a finite number of subsets $\Omega_j$ , $j=1,2,\ldots,J$ , so that: $\Omega=\bigcup_{j=1}^J\Omega_j,\qquad\Omega_j\cap\Omega_k=\varnothing\text{ if }j\ne k\tag{5.14}$ We let $\pi_j=P(x\in\Omega_j)$ and let $p_j(x)=p(x|x\in\Omega_j)$ , the probability distribution conditioned on $\Omega_j$ . We assume that the probabilities $\pi_j$ are known and that we can generate samples from the conditional distributions $p_j(x)$ . The sets $\Omega_j$ are called the strata. The partition theorem states:

Original page 24

$p(x)=\sum_{j=1}^J\pi_jp_j(x)\tag{5.15}$ $\mathcal{I}=\mathbb E[f(x)]=\sum_{j=1}^J\pi_j\mathbb E[f(x)|x\in\Omega_j]=\sum_{j=1}^J\pi_j\mathcal{I}_j\tag{5.16}$ if we let $\mathcal{I}_j=\mathbb E[f(x)|x\in\Omega_j]$ .

We now apply Monte Carlo to estimate each $\mathcal{I}_j$ , and the number of samples from each stratum need not be the same. Let $n_j$ be the number of samples from stratum $j$ , and we draw $n_j$ samples from the conditional distribution $p_j(x)$ : $x_1,\ldots,x_{n_j}\in\Omega_j$ . The Monte Carlo estimator for $\mathcal{I}_j$ is $\widehat{\mathcal{I}}_j=\frac1{n_j}\sum_{i=1}^{n_j}f(x_i)$ The overal Stratified Sampling (SS) estimator is $\widehat{\mathcal{I}}_{SS}=\sum_{j=1}^J\pi_j\widehat{\mathcal{I}}_j.\tag{5.17}$ It is easy to show that $\mathbb E[\widehat{\mathcal{I}}_{SS}]=\sum_{j=1}^J\pi_j\mathbb E[\widehat{\mathcal{I}}_j]=\mathbb E[f]\tag{5.18}$ namely $\widehat{\mathcal{I}}_{SS}$ is unbiased.

Variance of the estimator

Let: $\sigma_j^2=Var[f|\Omega_j]=\mathbb E[f^2|\Omega_j]-\mathcal{I}_j^2\tag{5.19}$ The quantity $\sigma_j^2$ is the variance of $X$ conditioned on $x\in\Omega_j$ .

The variance of $\widehat{\mathcal{I}}_{SS}$ is: $\operatorname{var}(\widehat{\mathcal{I}}_{SS})=\sum_{j=1}^J\frac{\pi_j^2\sigma_j^2}{n_j}\tag{5.20}$

Original page 25

which depends on the partition $\Omega_1,\ldots,\Omega_J$ .

Optimal allocation

A natural choice for $n_j$ is proportional allocation. Let $N$ be the total number of samples generated, then $n_j=\pi_jN$ . The variance becomes: $\operatorname{var}(\widehat{\mathcal{I}}_{SS})=\frac1N\sum_{j=1}^J\pi_j\sigma_j^2\tag{5.21}$ which is not the optimal choice of $n_j$ .

Alternatively we can determine $n_j$ by minimizing the variance. Namely we want to minimize the function: $F(n_1,n_2,\ldots,n_J)=\frac{a_1}{n_1}+\frac{a_2}{n_2}+\cdots+\frac{a_J}{n_J},\qquad a_j=\pi_j^2\sigma_j^2$ subject to the constraint: $n_1+n_2+\cdots+n_j=N.$ We solve it with the Lagrangian multiplier method. The Lagrangian $L$ is defined as: $L(x_1,x_2,\ldots,x_n,\lambda)=\frac{a_1}{x_1}+\frac{a_2}{x_2}+\cdots+\frac{a_n}{x_n}+\lambda(N-x_1-x_2-\cdots-x_n),$ where $\lambda$ is the Lagrange multiplier.

Take the partial derivatives of $L$ with respect to each $n_j$ and set them to zero: $\frac{\partial L}{\partial n_j}=-\frac{a_j}{n_j^2}-\lambda=0\quad\text{for }j=1,2,\ldots,J.$ Solve for $\lambda$ : $-\frac{a_i}{n_i^2}=\lambda\quad\Rightarrow\quad n_j^2=-\frac{a_i}\lambda\quad\Rightarrow\quad n_j=\sqrt{-\frac{a_j}\lambda}.$ From $n_1+n_2+\cdots+n_J=N$ we obtain: $\sum_{j=1}^J\sqrt{\frac{a_j}{-\lambda}}=N\quad=>\quad\lambda=-\left(\frac1N\sum_{j=1}^J\sqrt{a_j}\right)^2.$

Original page 26

Substitute $\lambda$ back into the expression for $n_j$ : $\begin{aligned} n_j&=\sqrt{-\frac{a_j}\lambda}=\sqrt{\frac{a_iN^2}{\left(\sum_{j=1}^J\sqrt{a_j}\right)^2}}\\ &=\frac{N\sqrt{a_j}}{\left(\sum_{j=1}^J\sqrt{a_j}\right)}=\frac{N\pi_j\sigma_j}{\left(\sum_{j=1}^J\pi_j\sigma_j\right)}. \end{aligned}$ This choice makes: $\operatorname{var}(\widehat{\mathcal{I}}_{SS})=\frac1N\left(\sum_{j=1}^J\pi_j\sigma_j\right)^2\tag{5.22}$

6 Sampling a univariate distribution

In the MC methods, a critical step is to draw samples from a given distribution. Starting from this section, we shall discuss how to do it. In this section, we consider the univariate distributions first.

6.1 Pseudorandom number generation

The foundation of random sampling methods is the ability to generate uniformly distributed random numbers, and many methods for sampling from other distributions are based on that. However, generating truly random numbers is very difficult in practice, and as such in most of the Monte Carlo simulations, pseudo-random numbers are usually used.

Pseudorandom numbers are numbers that mimic the behaviours of truly random numbers (in practise, this often means that they can pass a set of statistical tests). In other word, they appear to be random but are actually produced by a deterministic process. That deterministic process is called a pseudorandom number generator (PRNG).

The typical structure of a PRNG is as follows. There is a finite set $S$ of states, and a function $f:S\to S$ . There is an output space $U$ , and an output function $g:S\to U$ . We will always take the output space to be $[0,1]$ because we want to sample from $U[0,1]$ . The generator is given an initial value $S_0$ for the state, called the seed. The seed is typically provided

Original page 27

by the user. Then, a sequence of random numbers is generated by defining: $X_n=f(X_{n-1}),\qquad U_n=g(X_n),\qquad n=1,2,3,\ldots\tag{3.1}$ Note that there is nothing random here. If we generate a sequence of numbers with this procedure and then generate another sequence using the same seed, the two sequences will be identical. Since the state space is finite, $S_n$ must eventually return to a state it was in at some time earlier. The smallest integer $p$ such that for some state the function $f$ returns to that state after $p$ iterations is called the period of the generator. Obviously, longer periods are desirable. While noting that there are many PRNG available, we here present the simple linear congruential generators (LCG) as an example.

Linear Congruential Generator. It is defined by the recurrence relation: $X_{n+1}=(a\cdot X_n+c)\bmod m,\qquad U_n=X_n/m,$ where $a$ is the multiplier, $c$ is the increment, and $m$ is the modulus. The performance of LCG depends critically on the algorithm parameters $a,c$ and $m$ that should be carefully chosen.

Hereafter we do not distinguish the difference between the pseudo random numbers and the truly random ones. That is, we assume that we can generate random samples from a uniform distribution, $U[a,b]$ where $a<b$ are two constants.

6.2 The Box-Muller method

Next let us see how to sample a Gaussian distribution. The Gaussian distribution can be effectively sampled using a special algorithm – the Box-Muller method. Specifically, suppose that $u_1$ and $u_2$ are two independent samples drawn from the uniform distribution $U[0,1]$ . It follows that $x_1=\sqrt{-2\ln u_1}\sin(2\pi u_2),\qquad x_2=\sqrt{-2\ln u_1}\cos(2\pi u_2),\tag{6.1}$ follow standard Gaussian distribution and $x_1$ and $x_2$ are independent. This can be proved by considering the transform of random variables. Recall that, if we can write $u_1$ and $u_2$ as functions of $x_1$ and $x_2$ , i.e., $u_1=u_1(x_1,x_2)$ and $u_2=u_2(x_1,x_2)$ , the joint distribution of $x_1,x_2$ can be written as, $p(x_1,x_2)=p(u_1(x_1,x_2),u_2(x_1,x_2))|\det J(x_1,x_2)|,\tag{6.2}$

Original page 28

where $J$ is the Jacobian of the function $(u_1,u_2)$ . Because $p(u_1,u_2)=1$ (uniform distribution), it is easy to see that $p(x_1,x_2)=|detJ(x_1,x_2)|$ . Now we first derive from Eq. (6.1) $u_1=\exp(-x_1^2/2-x_2^2/2),\qquad u_2=\frac1{2\pi}\tan^{-1}(x_2/x_1),\tag{6.3}$ from which we can calculate the determinant of the Jacobian: $\begin{aligned} \det J&=\det\frac{\partial(u_1,u_2)}{\partial(x_1,x_2)}\\ &=\begin{vmatrix} \exp(-x_1^2/2-x_2^2/2)(-x_1)&\exp(-x_1^2/2-x_2^2/2)(-x_2)\\[4pt] -\frac1{2\pi}\frac{x_2}{x_1^2+x_2^2}&\frac1{2\pi}\frac{x_1}{x_1^2+x_2^2} \end{vmatrix}\\ &=-\frac1{\sqrt{2\pi}}\exp(-x_1^2/2)\frac1{\sqrt{2\pi}}\exp(-x_2^2/2), \end{aligned}\tag{6.4}$ and it follows directly that $p(x_1,x_2)=\frac1{\sqrt{2\pi}}\exp(-x_1^2/2)\frac1{\sqrt{2\pi}}\exp(-x_2^2/2)$ which implies that $x_1$ and $x_2$ each follow $N(0,1)$ and are independent of each other.

6.3 Inverse transform sampling

The Box-Muller method can only generate samples from Gaussian distributions. Now we consider sampling from an arbitrary univariate distribution. Namely let $x\in D\subset\mathbb R$ be a univariate random variable defined in a compact set $D$ and $F(x)$ be its cumulative density function. For simplicity we assume that $F$ is strictly increasing and strictly positive over $D$ . As a result function $F$ is invertible and its inverse is $F^{-1}:[0,1]\to D$ . The method of inverse sampling proceeds as follows,

first generate a set of samples $u_1,\ldots,u_n$ from $U[0,1]$ ;
for each $i=1\ldots n$ , let $x_i=F^{-1}(u_i)$ ;

and $x_1,\ldots,x_n$ are samples from the distribution $F(x)$ . To verify the correctness of this method, we need to show: let $u$ be a random variable following $U[0,1]$ , and let random variable $x=F^{-1}(u)$ , the CDF of $x$ is $F$ .

Original page 29

Proof: let $a$ be any number in $D$ , we have $P(x<a)=P(F^{-1}(u)<a)=P(u<F(a))=F(a)$ as $u$ follows the uniform distribution.

Note: a graph will be shown here on how to draw samples.

We note here that though conceptually simple, the implementation of the method can be quite expensive, as computing the inversion of an arbitrary CDF can be quite difficult and sometimes requires solving a nonlinear equation for each sample.

6.4 Rejection sampling

The rejection sampling method applies to multivariate distributions in principle. However, because its performance is severely limited by dimensionality, it is generally not preferred for high dimensional problems. In this section we restrict ourselves in the univariate case. The rejection sampling method relies on the availability of a proposal distribution. Namely suppose we want to generate samples from a target distribution $p$ , and we choose a proposal distribution $q$ from which we can sample. Now we choose a constant $c$ such that $c\geq\sup_x\frac{p(x)}{q(x)}$ and the rejection sampling method generates a sample $x$ from the target distribution $p$ using the following procedure,

generate a number $u$ from $U[0,1]$ ;
generate a number $x'$ from $q$ ;
let $x=x'$ if $u<\frac{p(x)}{cq(x)}.$

Next we show that the samples drawn this way follow the target distribution. Specifically we look at the cumulative density function, for any $x^*$ in the

Original page 30

state space, $\begin{aligned} P(x<x^*)&=P(x'<x^*|x'\text{ accepted})=P\left(x'<x^*\middle|u<\frac{p(x')}{cq(x')}\right)\\ &=\frac{P\left(x'<x^*,u<\frac{p(x')}{cq(x')}\right)}{P\left(u<\frac{p(x')}{cq(x')}\right)} \end{aligned}\tag{6.5}$ For the numerator, using the feature of conditional probability, we play a trick of introducing an artificial random variable $x''$ which follows $q(x'')$ : $\begin{aligned} &P\left(x'<x^*,u<\frac{p(x')}{cq(x')}\right)\\ &=\int P\left(x'<x^*,u<\frac{p(x')}{cq(x')}\middle|x'=x''\right)q(x'')dx''\\ &=\int P\left(x''<x^*,u<\frac{p(x'')}{cq(x'')}\right)q(x'')dx''\\ &=\int_{\{x''<x^*\}}P\left(u<\frac{p(x'')}{cq(x'')}\right)q(x'')dx''\\ &=\int_{\{x''<x^*\}}\frac{p(x'')}{cq(x'')}q(x'')dx''=\frac1c\int_{-\infty}^{x^*}p(x'')dx'' \end{aligned}\tag{6.6}$ Note that in the 4-th equality in the equation above we use the fact that $c\geq\sup\frac{p(x'')}{q(x'')}$ to ensure that $\frac{p(x'')}{cq(x'')}\leq1$ . For the denominator we use the similar trick $\begin{aligned} P\left(u<\frac{p(x')}{cq(x')}\right)&=\int P\left(u<\frac{p(x')}{cq(x')}\middle|x'=x''\right)q(x'')dx''\\ &=\int P\left(u<\frac{p(x'')}{cq(x'')}\right)q(x'')dx''\\ &=\int\frac{p(x'')}{cq(x'')}q(x'')dx''=\frac1c\int p(x'')dx''=\frac1c \end{aligned}\tag{6.7}$ Substituting Eqs (6.6) and (6.7) into Eq. (6.5), we obtain, $P(x<x^*)=\int_{-\infty}^{x^*}p(x'')dx'',$ which implies that $x$ follows the target distribution $p(\cdot)$ . Finally we need comment on the choice of the constant $c$ . As is mentioned earlier, for the method to work, we need to have $c\geq\sup_x\frac{p(x)}{q(x)}$ . On the other hand, from an efficiency perspective, a larger $c$ means a higher chance that a proposed sample gets rejected¹ which reduces the sampling efficiency. Therefore ideally an optimal value of $c$ should be taken to be $c=\sup_x\frac{p(x)}{q(x)}$ .

Original page 31

7 Markov Chain Monte Carlo

7.1 Basics of Markov Chains

Starting from this section, we will consider the problem of generating samples from an arbitrary distribution especially for multivariate distributions. To this end, one of the most popular methods is Markov chain Monte Carlo (MCMC). Before discussing MCMC, in this section we will review the Markov Chain Process, which is the foundation of MCMC. For simplicity we will only discuss the Markov Chains in a finite state space, which means that the stochastic process only takes values from a finite set, while noting that the properties discussed here apply to those defined in a continuous state space.

Markov chain. Let $x_t$ be a stochastic process defined on a finite set $\{x_1,\ldots,x_s\}$ , and $x_t$ is called a Markov process if $p(x_t|x_{t-1},\ldots,x_1)=p(x_t|x_{t-1}).$ Loosely speaking, for a Markov chain process, the distribution of future outcomes only depends on its present state. Suppose that the distribution of $x_{t-1}$ is $\pi_{t-1}(\cdot)$ , and that of $x_t$ is $\pi_t(\cdot)$ , it is easy to see that $\sum_{x_{t-1}\in S}\pi_{t-1}(x_{t-1})p(x_t|x_{t-1})=\pi_t(x_t).\tag{7.1}$ Homogeneous. The chain is called homogeneous if $p(x_t|x_{t-1})=p(x_{t'}|x_{t'-1})$ for any $t$ and $t'$ , with $\sum_{x_t\in S}p(x_t|x_{t-1})=1$ for any $t$ .

In practice, the conditional distribution can be represented as a transition matrix $T^t=\{p_{i,j}\}_{i,j=1}^s$ where $p_{i,j}$ represents the probability that the chain jumps from state $x_i$ to the state $x_j$ . In the homogeneous case, $T^t=T$ . Eq. (7.1) can be writen as $\pi_{t-1}T=\pi_t$ .

Original page 32

Example: suppose the state space is $\{x_1,x_2,x_3\}$ and a transition matrix is $T=\begin{pmatrix}0&1&0\\0&0.1&0.9\\0.6&0.4&0\end{pmatrix},$ and note that the elements in each row add up to 1.

Invariant distribution. Some homogeneous Markov chain processes admit an invariant distribution (or stationary distribution or equilibrium distribution) which is defined as a distribution $\pi$ satisfying $\pi T=\pi$ . That is, if $x_t$ follows distribution $\pi$ , so does $x_{t+1}$ . A very important property of a Markov chain is that the chain admits an invariant distribution and converges to its invariant distribution as $t$ increases, if the chain satisfies the detailed balance condition: $\pi(x_t)T(x_{t-1}|x_t)=\pi(x_{t-1})T(x_t|x_{t-1}).\tag{7.2}$ Summing both sides of Eq. (7.2) over $x_{t-1}\in S$ yields $\pi(x_t)=\sum_{x_{t-1}\in S}\pi(x_{t-1})T(x_t|x_{t-1}),$ namely $\pi=\pi T$ .

7.2 Metropolis-Hastings

Recall that our goal is to draw samples from a given distribution. Now suppose the given distribution is $\pi$ , and the main idea of MCMC is to design a Markov chain which has the target distribution $\pi$ as its invariant distribution. That means, if we draw the chain sufficiently long, the distribution of $x_t$ will converge to the target distribution, and we can regard $x_t$ as a sample drawn from the target distribution. In practice MCMC is usually implemented via the Metropolis-Hasting (MH) scheme.

As is mentioned earlier, Metropolis-Hasting algorithm (MH) provides a way to construct such a Markov chain. More precisely, the MH algorithm constructs the chain using the following procedure,

initial at $x_0$
For $t=0$ to $n-1$

Original page 33

draw $x^*$ from a proposal distribution $q(x^*|x_t)$ ;
$a=\min\{1,\frac{\pi(x^*)}{\pi(x_t)}\frac{q(x_t|x^*)}{q(x^*|x_t)}\}$
draw $u\sim U[0,1]$ ;
if $u<a$ , let $x_{t+1}=x^*$ , otherwise let $x_{t+1}=x_t$ .

In the MH procedure, $a$ is often called the acceptance probability as that the algorithm essentially accepts the new proposed point $x^*$ with probability $a$ . Here we omit the proof and just state that the Markov chain constructed using the MH procedure has the target distribution $\pi(\cdot)$ as its invariant distribution. As one can see here a key in the MH algorithm is the choice of the proposal distribution $q(\cdot|\cdot)$ . Two popular choices of the proposal distributions are random walks around $x_{t-1}$ : $q(x^*|x_t)=N(x_t,\sigma^2I),$ and $q(x^*|x_t)=U[x_t-\sigma,x_t+\sigma]^d,$ where in both proposals the parameter $\sigma$ plays the role as the stepsize. It should be clear that the first proposal is a Gaussian random walk centered at $x_t$ . In the second proposal $U[x_t-\sigma,x_t+\sigma]^d$ means a uniform distribution defined on the $d$ -dimensional hypercube $[x_t-\sigma,x_t+\sigma]$ , and so it is a uniform random walk around $x_t$ . It should be noted here that these two proposals are “symmetric” in the sense that $q(x'|x)=q(x|x').$ and as a result, when either of them is used in the MH algorithm, the acceptance probability becomes $a=\min\left\{1,\frac{\pi(x^*)}{\pi(x_{t-1})}\right\}.$ It is also important to mention here that, a remarkable feature of MCMC is that the method does not require the exact knowledge of the target distribution; rather one can draw samples from a distribution if it is known up to a constant. That is, if a distribution is in the form of $\pi(x)=cg(x)$ where $g(x)$ is known but $c$ is an unknown constant, one can still use MCMC to draw samples from it (why?).

Example and excercise: say the target distribution is standard Gaussian: $\pi(x)=\frac1{\sqrt{2\pi}}\exp(-x^2/2)$ ; suppose that we only know the part $\exp(-x^2/2)$

Original page 34

and do not know the constant $1/\sqrt{2\pi}$ . Try to write a MH algorithm to generate samples from it using the uniform random walk proposal with $\sigma=0.5$ .

Exercise: target distribution $\pi(x)=c\exp(-(1-x)^2-10(y-x^2)^2)$ where $c$ is an unknown constant. Write a MH algorithm to generate samples from the target distribution using a Gaussian random walk proposal. Show the scatter plots of the points generated. Try different values of $\sigma$ (the stepsize) and compare the results.

Remark: choosing the stepsize is also a critical issue in MH. If the stepsize is too small, the chain is moving very slow, and if the stepsize is too large the new point is very likely to be rejected. On the other hand, in this case we can also see that, for MH with a random walk proposal, very high acceptance probability is actually not desired (obviously very low acceptance probability is not desired as well), as that often means the stepsize is too small.

7.3 Independent sampler

A special implementation of the MH algorithm is the so-called independent sampler. In this method, the proposal distribution is taken to be independent of $x_t$ , namely $q(x^*|x_t)=q(x^*)$ , and as a result the acceptance probability becomes, $a=\min\left\{1,\frac{\pi(x^*)}{\pi(x_t)}\frac{q(x_t)}{q(x^*)}\right\}.$ One can see here that if $q(\cdot)\approx\pi(\cdot)$ , the acceptance probability is close to 1; that is all the proposed points will be accepted. As such the independent sampler is usually used when one can construct an approximation of the target distribution, which is used as the proposal.

It should be noted here that, unlike the random walk MH, in the independent sampler, the higher the acceptance probability is, the better the algorithm performs, as that often means the proposal is very close to the target distribution.

Original page 35

7.4 Gibbs sampling

The Gibbs sampler is a special implementation of the MH algorithm which is used in a special class of problem. Namely, suppose $x$ is a $d$ -dimensional random variable and is written as $x=(x_1,\ldots,x_d),$ and follows distribution $\pi(x_1,\ldots,x_d)$ . We introduce the notation $x_{-j}=(x_1,\ldots,x_{j-1},x_{j+1},\ldots,d).$ The Gibbs sampler can be used if we can directly draw samples from the conditional distributions $\pi(x_j|x_{-j})$ for $i=1\ldots d$ . Specifically at step $t$ , suppose that the state obtained from the previous step is $x^t$ , and the Gibbs sampler uses the following proposal: for any $j=1\ldots d$ , $q(x^*|x^t)=\begin{cases}\pi(x_j^*|x_{-j}^t),&\text{if }x_{-j}^*=x_{-j}^t\\0,&\text{otherwise}.\end{cases}\tag{7.3}$ The proposal can be understood as the following: suppose the current state is $x^t=(x_1^t,\ldots,x_d^t)$ , the new sample drawn from the proposal $q(x^*|x^t)$ is $x^*=(x_1^t,\ldots,x_{j-1}^t,x_j^*,x_{j-1}^t,\ldots,x_d^t),$ where $x_j^*$ is drawn from $\pi(x_j^*|x_{-j}^t)$ . Remarkably we can prove that the acceptance probability with such a proposal is 1: $\begin{aligned} a(x^*,x^t)&=\min\left\{1,\frac{\pi(x^*)}{\pi(x^t)}\frac{q(x^t|x^*)}{q(x^*|x^t)}\right\}\\ &=\min\left\{1,\frac{\pi(x_1^t,\ldots,x_{j-1}^t,x_j^*,x_{j+1}^t,x_d^t)}{\pi(x_1^t,\ldots,x_{j-1}^t,x_j^t,x_{j+1}^t,x_d^t)}\frac{\pi(x_j^t|x_{-j}^*)}{\pi(x_j^*|x_{-j}^t)}\right.\\ &=\min\left\{1,\frac{\pi(x_j^*|x_{-j}^t)\pi(x_{-j}^t)}{\pi(x_j^t|x_{-j}^t)\pi(x_{-j}^t)}\frac{\pi(x_j^t|x_{-j}^t)}{\pi(x_j^*|x_{-j}^t)}\right\}\\ &=1, \end{aligned}$ which means the samples drawn from this proposal will all be accepted.

It is important to note here, in principle, in each step one can choose any $1\leq j\leq d$ and use the Gibbs sampler to draw a new sample (this new sample only differs from the current one at the $j$ -th compopent). However for effectiveness, we usually chose to generate a new sample by updating all the components from $i=1\ldots d$ , which results in the following algorithm of the Gibbs sampler:

Original page 36

choose a $x^0$ ;
for $t=0:(T-1)$ do
- draw $x_1^{t+1}$ from $\pi(x_1|x_{-1}^t)$
- ......
- draw $x_j^{t+1}$ from $\pi(x_j|x_{-j}^t)$
- ......
- draw $x_d^{t+1}$ from $\pi(x_d|x_{-d}^t)$ ,

and in practice a variate of the algorithm is often used:

choose a $x^0$ ;
for $t=0:(T-1)$ do
- let $x^{t+1}=x^t$ ;
- draw $x_1^{t+1}$ from $\pi(x_1|x_{-1}^t)$ ;
- ......
- draw $x_j^{t+1}$ from $\pi(x_j|x_{-j}^t)$ ;
- ......
- draw $x_d^{t+1}$ from $\pi(x_d|x_{-d}^t)$ .

Question: what is the difference between these two algorithms?

7.5 Adaptive Metropolis

This section is optional. Recall that, in the standard MH algorithm the proposal is often taken to be “isotropic” Gaussian, $q(x^*|x_t)=N(x_t,\sigma^2I)$ . However, it has been known that it is often more efficient to using the following proposal $q(x^*|x_t)=N(x_t,\sigma^2\Sigma)$ where $\Sigma$ is the covariance of the target distribution. It should be clear that direct use of this proposal is not possible because we do not know the target covariance $\Sigma$ in advance. One idea is to adaptively estimate the covariance matrix from the sample history, and the adaptive Metropolis algorithm is

Original page 37

based on this idea. The basic scheme of AM is quite straightforward: suppose at step $t$ , we have already generated a set of samples $x_1,\ldots,x_t$ , we first estimate the sample covariance as $\hat\Sigma_t=\frac1{t-1}\sum_{i=1}^t(x_i-\hat x_t)(x_i-\hat x_t)^T,\qquad\hat x_t=\frac1t\sum_{i=1}^tx_i,\tag{7.4}$ where the superscript $T$ means matrix transpose, and update the proposal for the next step as $q(x^*|x_t)=N(x_t,\sigma^2\hat\Sigma_t).$ Recursive formula. It should be noted here that computing the sample covariance using Eq. (7.4) can be rather time consuming especially when the number of samples is large, and thus in the numerical implementation, we use a recursive formula: namely suppose that for a set of samples $\{x_1,\ldots,x_{t-1}\}$ , we have already obtained its sample mean and sample variance which are respectively $\hat x_{t-1}$ and $\Sigma_{t-1}$ , and as a new sample $x_t$ becomes available, the sample mean and variance of $\{x_1,\ldots,x_{t-1},x_t\}$ can be calculated as $\hat x_t=\frac{t-1}t\hat x_{t-1}+\frac1tx_t,$ $\hat\Sigma_t=\frac{t-2}{t-1}\hat\Sigma_{t-1}+\frac1t(x_t-\hat x_{t-1})(x_t-\hat x_{t-1})^T.$ The complete algorithm of AM for drawing samples from a target distribution $\pi(x)$ proceeds as follows:

initialization: choose a starting point $x_0$ , a stepsize $\sigma$ and an initial covariance matrix $\hat\Sigma_0$ ;
standard MH: generate $n_0$ samples $x_1,\ldots,x_{n_0}$ from the target distribution $\pi(x)$ using proposal $q(x^*|x_t)=N(x_t,\sigma^2\Sigma_0)$ ;
estimate the sample mean and variance of $\{x_1,\ldots,x_{n_0}\}$ : $\begin{gathered} \hat\Sigma_{n_0}=\frac1{n_0-1}\sum_{i=1}^t(x_i-\hat x_{n_0})(x_{n_0}-\hat x_{n_0})^T,\\ \hat x_{n_0}=\frac1{n_0}\sum_{i=1}^{n_0}x_i, \end{gathered}$
for $t=n_0$ to $n-1$ do
- draw $x^*$ from the proposal $q(x^*|x_t)=N(x_t,\sigma^2\Sigma_t)$ ;
- $a=\min\{1,\frac{\pi(x^*)}{\pi(x_{t-1})}\frac{q(x_t|x^*)}{q(x^*|x_t)}\}$
- draw $u\sim U[0,1]$ ;

Original page 38

- if $u<a$ , let $x_{t+1}=x^*$ , otherwise let $x_{t+1}=x_t$ ;
- update the sample mean and variance: $\begin{gathered} \hat x_t=\frac{t-1}t\hat x_{t-1}+\frac1tx_t,\\ \hat\Sigma_t=\frac{t-2}{t-1}\hat\Sigma_{t-1}+\frac1t(x_t-\hat x_{t-1})(x_t-\hat x_{t-1})^T. \end{gathered}$

Finally it should be noted that the AM algorithm is essentially non-Markovian as its present proposal $q(x^*|x_t)$ depends on the whole history $\{x_1,\ldots,x_t\}$ . Nevertheless it can still be proved that the chain generated by AM converges to the target distribution. The proof is different from that for the standard MH algorithm which is based on the detailed balance condition, and will not be discussed here.

7.6 Diagnostics of MCMC algorithms

This section is optional. In this section we discuss the tools for evaluating the performance of MCMC algorithms. A major limitation of MCMC is that the samples generated are correlated, and on the other hand, in MC simulation we usually prefer to use independent samples. To this end, it is desired that the series correlation of the chain should decaying as fast as possible and in practice, we often resort to the auto-correlation function (ACF) of the samples to evaluate the performance of the method. Mathematically the ACF of a stationary stochastic process (where most MCMC sampels are) $x_t$ is defined as, $\operatorname{ACF}(\tau)=\frac{\mathop{\mathrm{Cov}}[x_t,x_{t+\tau}]}{\mathop{\mathrm{Var}}[x_t]},$ where $\tau$ is the time lag, i.e. the time difference between the two points. We note here that in many popular programming languages, sample ACF can be calculated directly (for example, in Matlab it can be computed by using the function $autocorr$ ). The ACF plots show how the sample autocorrelation between the terms of the chain decreases as a function of their lag and loosely speaking we prefer that the ACF decays fast to zero. Using the ACF, we can compute another popular performance measure of MCMC, the effective sample size (ESS). Roughly speaking ESS measures effectively how many effectively independent samples that a set of MCMC samples are equivalent to. ESS can be calculated with the following empirical formula: $N_{eff}=\frac{N}{1+2\sum_{\tau=0}^\infty\operatorname{ACF}(\tau)},$ where $N$ is the number of MCMC samples.

Original page 39

8 Nonlinear linear regression and optimization

8.1 Linear and nonlinear regression

Linear Regression. Linear regression is a statistical method used to model the relationship between a dependent variable $y$ and some independent variables $x=(x_1,\ldots,x_d)$ . It assumes a linear relationship between the variables. Specifically the relationship between $y$ and $x$ can be expressed as: $y=\theta_0+\theta_1x_1+\theta_2x_2+\cdots+\theta_dx_d+\epsilon$ where $\theta_0,\theta_1,\ldots,\theta_d$ are the coefficients and $\epsilon$ is the error term.

Suppose that we have a set of $n$ ( $n>d$ ) paired data points $\{(x_i,y_i)\}_{i=1}^n$ , and we will compute the coefficients $\theta_0,\ldots,\theta_d$ by minimizing the sum of squared residuals: $\operatorname{RSS}=\sum_{i=1}^n(y_i-(\theta_0+\theta_1x_1^i+...+\theta_dx_d^i))^2$ This is often called Ordinary Least Squares (OLS) and the solution can be found analytically. While the linear regression is computationally simple, it assumes a linear relationship, which may not hold for all datasets.

Nonlinear Regression. Nonlinear regression models data where the relationship between the dependent and independent variables is not linear. The general form of a nonlinear model is: $y=f(x,\theta)+\epsilon$ where $f$ is a nonlinear function of $x$ and $\theta$ (parameters). Here $\theta=(\theta_0,\theta_1,\ldots,\theta_p)$ and note that in generally $p\ne d$ .

Nonlinear Least Squares (NLS): The parameters ( $\theta$ ) are estimated by minimizing the Residual Sum of Squares (RSS): $\operatorname{RSS}=\sum_{i=1}^n[y_i-f(x_i,\theta)]^2.$ Since the equations are nonlinear in $\theta$ , iterative optimization methods are used.

Original page 40

Examples of Nonlinear Models

Exponential Model: $y=\theta_0e^{\theta_1x}+\epsilon$
Logistic Model: $y=\frac{\theta_0}{1+e^{-(\theta_1+\theta_2x)}}+\epsilon$
Power-law Model: $y=\theta_0x^{\theta_1}+\epsilon$
Neural Network Model: $y=a(Wx+b)+\epsilon,$ where $a$ is a nonlinear activation function and a common example is the Sigmoid Activation: $a(z)=1/(1+e^{-x})$ .

8.2 Basics of Optimization Algorithms

In this section we give a very brief overview of some basic algorithms for unconstrained optimization. An unconstrained optimization is simply $\min_\theta F(\theta)$ for a given real-valued function $F(\theta)$ . There are many (iterative) algorithms can be used to solve the unconstrained optimization problem. One of the most common algorithms is the gradient descent method, which conducts the following iteration in each step: $\theta_{k+1}=\theta_k-\alpha_k\nabla F(\theta_k)$ where $\alpha_k$ is the step size and $\nabla F(\theta_k)$ is the gradient of $F$ at point $\theta_k$ . The gradient descent method is known to converge linearly to the optimal solution under certain conditions. Loosely speaking linear convergence means that $\frac{\|\theta_{k+1}-\theta^*\|}{\|\theta_k-\theta^*\|}\leq C$ where $x^*$ is the optimal solution and $C$ is a constant in $(0,1)$ . If high-order methods are desired, the Newton’s method is often used: $\theta_{k+1}=\theta_k-\alpha_kH^{-1}(\theta_k)\nabla F(\theta_k)$

Original page 41

where $H(\cdot)$ is the Hessian matrix of $F$ : $H=\begin{bmatrix} \frac{\partial^2F}{\partial\theta_1^2}&\frac{\partial^2F}{\partial\theta_1\partial\theta_2}&\cdots&\frac{\partial^2F}{\partial\theta_1\partial\theta_p}\\[5pt] \frac{\partial^2F}{\partial\theta_2\partial\theta_1}&\frac{\partial^2F}{\partial\theta_2^2}&\cdots&\frac{\partial^2F}{\partial\theta_2\partial\theta_p}\\ \vdots&\vdots&\ddots&\vdots\\[3pt] \frac{\partial^2F}{\partial\theta_p\partial\theta_1}&\frac{\partial^2F}{\partial\theta_p\partial\theta_2}&\cdots&\frac{\partial^2F}{\partial\theta_p^2} \end{bmatrix}.$ It is well-known that the Newton’s method is (under certain conditions) quadratically convergent to the optimal solution, in the sense that $\frac{\|\theta_{k+1}-\theta^*\|}{\|\theta_k-\theta^*\|^2}\leq C$ where $C$ is a constant in $(0,1)$ .

8.3 Optimization algorithms for NLS

Here we start by applying optimization algorithms method to the NLS problem introduced in Section 8.1. $\begin{gathered} F(\theta)=\sum_{i=1}^n[y_i-f(x_i,\theta)]^2=(r_i(\theta))^2,\\ r_i(\theta)=y_i-f(x_i,\theta). \end{gathered}$ If we want to use the gradient descent method, we need the gradient of the function $F$ . Now we introduce the notation: $r(\theta)=\begin{bmatrix}r_1(\theta)\\r_2(\theta)\\\vdots\\r_n(\theta)\end{bmatrix}$ and we can write $F(\theta)$ as $F(\theta)=r(\theta)^Tr(\theta).$ We can compute the gradient of $F$ $\nabla F=\begin{bmatrix}\frac{\partial F}{\partial x_1}(\theta)\\[4pt]\frac{\partial F}{\partial\theta_2}(\theta)\\\vdots\\[3pt]\frac{\partial F}{\partial\theta_p}(\theta)\end{bmatrix}$

Original page 42

where $\frac{\partial F}{\partial\theta_j}(\theta)=\sum_{i=1}^n\frac{\partial r_i}{\partial\theta_j}(\theta)r_i(\theta),$ and also note that $\frac{\partial r_i}{\partial\theta_j}(\theta)=-\frac{\partial f(x_i,\theta)}{\partial\theta_j}(\theta).$ Because $r(\theta)$ is vector-valued function, we can compute its Jacobian as $J(\theta)=\begin{pmatrix} \frac{\partial r_1}{\partial\theta_1}(\theta)&\cdots&\frac{\partial r_1}{\partial\theta_p}(\theta)\\ ...&...&...\\ \frac{\partial r_n}{\partial\theta_1}(\theta)&\cdots&\frac{\partial r_n}{\partial\theta_p}(\theta) \end{pmatrix}$ It can be verified that $\nabla f(\theta)=J(\theta)^Tr(\theta)$ .

Gradient descent method: $\theta_{k+1}=\theta_k-\alpha_kJ(\theta_k)^Tr(\theta_k)$ In many cases, Newton’s method is preferred thanks to its quadratic convergence. If we want to use the Newton’s method, we need both the gradient and the Hessian matrix of $F(\theta)$ . Hessian matrix $H(\theta)$ is symmetric: $H=\begin{bmatrix} \frac{\partial^2F}{\partial\theta_1^2}&\frac{\partial^2F}{\partial\theta_1\partial\theta_2}&\cdots&\frac{\partial^2F}{\partial\theta_1\partial\theta_p}\\[5pt] \frac{\partial^2F}{\partial\theta_2\partial\theta_1}&\frac{\partial^2F}{\partial\theta_2^2}&\cdots&\frac{\partial^2F}{\partial\theta_2\partial\theta_p}\\ \vdots&\vdots&\ddots&\vdots\\[3pt] \frac{\partial^2F}{\partial\theta_p\partial\theta_1}&\frac{\partial^2F}{\partial\theta_p\partial\theta_2}&\cdots&\frac{\partial^2F}{\partial\theta_p^2} \end{bmatrix}.$ and its $(j,j')$ element is $\frac{\partial^2F}{\partial\theta_j\partial\theta_{j'}}= \sum_{i=1}^n\frac{\partial r_i}{\partial\theta_j}(\theta)\frac{\partial r_i}{\partial\theta_{j'}}(\theta) +\sum_{i=1}^nr_i(\theta)\frac{\partial^2r_i}{\partial\theta_j\partial\theta_{j'}}(\theta)$ Thus $H(\theta)=J(\theta)^TJ(\theta)+S(\theta)$ where $S(\theta)$ is formed by $S(j,j')=\sum_{i=1}^nr_i(\theta)\frac{\partial^2r_i}{\partial\theta_j\partial\theta_{j'}}(\theta).$

Original page 43

Newton’s method: $\theta_{k+1}=\theta_k-(J(\theta_k)^TJ(\theta_k)+S(\theta_k))^{-1}J(\theta_k)^Tr(\theta_k)$ Newton’s method involves the calculation of second order derivative of $r(\theta)$ , which in turn involves the second order derivative of the nonlinear regression model $f(x,\theta)$ with respect to $\theta$ . In practice, when the regression model is complex, its Hessian can be very difficult to compute. In this case, a simplified version of the Newton’s method is often used. It is easy to see that the term in the Hessian that involves the 2nd order derivative of $r(\theta)$ is $S(\theta)$ . As such we simply remove the term $S(\theta)$ in the Newton’s method, yielding the so-called Gauss-Newton method.

Gauss-Newton’s method: $\theta_{k+1}=\theta_k-(J(\theta_k)^TJ(\theta_k))^{-1}J(\theta_k)^Tr(\theta_k)$ Remark 1: In practice, one often does not compute the inverse matrix $(J(\theta_k)^TJ(\theta_k))^{-1}$ . Instead, one computes the update by solving a linear equation system: $\theta_{k+1}=\theta_k-\Delta\theta,$ with $(J(\theta_k)^TJ(\theta_k))\Delta\theta=J(\theta_k)^Tr(\theta_k).$ Remark 2: A potential problem is that $J(x)^TJ(x)$ is not invertible. A solution in this case is to add a small diagonal matrix $\delta I_d$ to $J(x)^TJ(x)$ , making it invertible: $\theta_{k+1}=\theta_k-(J(\theta_k)^TJ(\theta_k)+\delta I_d)^{-1}J(\theta_k)^Tr(\theta_k)$ $\theta_{k+1}=\theta_k-\Delta\theta,$ with $(J(\theta_k)^TJ(\theta_k)+\delta I_d)\Delta\theta=J(\theta_k)^Tr(\theta_k).$ Here $\delta$ is a small positive number and $I_d$ is the $d\times d$ identity matrix. This method is called the Levenberg-Marquardt method.