Topics for Nov 11 and 13

36-711 Fall 2002

Outline

1. Numerical integration and standard Bayesian calculations.

(a) Numerical integration to compute posterior probabilities, e.g. CI’s.  bayes.s

2. Sketch of some standard univariate numerical integration methods.

3. Sampling to produce pictures of posterior densities

(a) Sampling generalities.

(b) Exploring posterior distributions in Bayesian inference, using Monte Carlo. sampling.s

[later we will talk briefly about multivariate numerical integration]

1 Numerical integration and standard Bayesian calculations

Rarely will the anti-derivatives composed of elementary functions like x^{k}, e^{x}, sin(x), etc., that we learned in Calculus, be sufficient to compute definite integrals in real life.

We have already seen in the E-M algorithm that the integral

Q(θ, θ^{k}) = ∫ log f(X|θ)f(X|Y, θ^{k})dX

is a basic ingredient in maximum likelihood problems with missing data. In “simple” problems this can be evaluated in closed form, but usually it cannot and must somehow be approximated.

Another motivation for thinking about numerical integration methods is is in calculations for Bayesian inference. Mathematically similar integration problems arise in statistical inference for mixture models.

Given a (possibly multidimensional) parameter θ with prior distribution π(θ) and likelihood L(y|θ) = ∏_{i} L_{i}(y_{i}|θ), Bayesian inference is based on functionals (i.e. numerical functions) of the posterior distribution

p(θ|y) = \frac{L(y|θ)π(θ)}{∫ L(y|θ)π(θ)dθ}

Common functionals are

• The posterior probability of an interval [a, b], ∫ 1_{{θ∈[a,b]}}p(θ|y)dθ; • The posterior mean E[Θ|y] = ∫ θp(θ|y)dθ;

• The posterior variance Var(Θ|y) = ∫ (θ − E[Θ|y])^{2}p(θ|y)dθ;

• The “entire” posterior density p(θ|y).

A direct attack on the first three quantities requires at the very least two integrals:

E[g(θ)|y] = ∫ g(θ)p(θ|y)dθ = \frac{∫ g(θ)L(y|θ)π(θ)dθ}{∫ L(y|θ)π(θ)dθ}  (1)

The first part of our discussion will treat some ways of accomplishing the direct attack.

Paradoxically, we can often approximate p(θ|y) quite well without evaluating any integrals, by doing Monte Carlo simulation: we generate iid random draws from p(θ|y), and then draw a histogram or a kernel density estimate. This also permits an end-run around the integral calculations for the first three problems; for, by the Law of Large Numbers

∫ g(θ)p(θ|y)dθ ≈ \frac{1}{n} ∑_{1}^{n} g(θ_{i})

where θ_{1}, …, θ_{n} is a random sample from p(θ|y).

Often it is difficult to sample directly from p(θ|y). The second part of our discussion will focus on elementary ways to indirectly generate a random sample from p(θ|y), by first obtaining a sample from π(θ) and then adjusting it somehow.

2 Sketches of some standard numerical integration methods

In addition to Monahan’s text, the text Elements of Statistical Computing by Thisted is a good reference here.

2.1 Newton-Cotes rules

The simple methods one learns about in Calculus can actually be quite useful. We begin with a “nice” function f(x) on a fixed interval (a, b) and an increment h = (b−a)/n, so that a = a+0h, …, x_{i} = a+ih, …, b = a + nh.

• Rectangle rule.

∫_{a}^{b} f(x)dx ≈ h ∑_{i=1}^{n−1} f(x_{i})

– Open vs Closed: Are the endpoints included in the calculation (maybex_{i} = a + (n + \frac{1}{2})?). Most of the rules below have “open” and “closed” versions. We will mainly discuss the “closed” ones, since they are a little simpler to describe.

– Extended/composite rule: Add the answers

hf(x_{i}),

which are exact as long as f(x) is constant on each subinterval.

• Trapezoid rule.

∫_{a}^{b} f(x)dx ≈ \frac{h}{2}f(a) + h ∑_{i=1}^{n−1} f(x_{i}) + \frac{h}{2}f(b)

– Can show the approximation is exact up to O(f′′(x^{∗})/n^{2}), so doubling the number of points quadruples the accuracy of the integral approximation.

– Extended/composite rule: Add the answers

\frac{h}{2} [f(x_{i}) + f(x_{i+1})]

which are exact as long as f(x) is piecewise linear on the intervals defined by x_{0}, …, x_{n}.

• A reasonably efficient trapezoid rule.

trap.int _function(fcn, a, b, maxcalls, releps) {
#
# Trapezoid rule, doubling number of points in the interval until
# estimated integral stops changing values.  Re-uses old function
# calls, so that no redundant calls to the integrand (fcn) are made.
#
#
k _ceiling(log(maxcalls - 1)/log(2))
h _b-a
n _1
I0 _h*(fcn(b) - fcn(a))/2
I1 _(I0+1)/(releps + 1)
j _1
while ((j<=k) && (abs(I1-I0)>releps*abs(I0))) {
         tmp _a + h/2
         p _0
         for (i in (0:(n-1))) {
                  p _p + fcn(tmp + i*h)
         }
         I0 _I1
         I1 _(I0+h*p)/2
         h _h/2
         n _2*n
         j _j + 1
}
return(I1,j,k)
}

• Simpson’s Rule.

∫_{a}^{b} f(x)dx ≈ \frac{h}{3} [f(a) + 4f(x_{1}) + 2f(x_{2}) + 4f(x_{3}) + 2f(x_{4}) + ⋯ + f(b)]

– Can show the approximation is exact up to O(f^{(4)}(x^{∗})/n^{4}): doubling the number of points increases the the accuracy of the integral approximation by a factor of 16.

– Extended/composite rule: Add the answers

\frac{h}{2} [f(x_{i}) + 4f(x_{i+1}) + f(x_{i+2})]

which are exact as long as f(x) is piecewise linear on the each pair of adjacent subintervals defined by x_{0}, …, x_{n}.

• General Newton-Cotes methods.

– They evaluate the function on a grid of equally-spaced points.

– They give the exact answer when integrating certain polynomials.

Clearly this is just a linear system to solve:

w_{1}f^{0}(x_{1}) + w_{2}f^{0}(x_{2}) + ⋯ + w_{k}f^{0}(x_{k}) = I_{0} ; w_{1}f^{1}(x_{1}) + w_{2}f^{1}(x_{2}) + ⋯ + w_{k}f^{1}(x_{k}) = I_{1} ; w_{1}f^{2}(x_{1}) + w_{2}f^{2}(x_{2}) + ⋯ + w_{k}f^{2}(x_{k}) = I_{2} ; ⋮ ; w_{1}f^{k−1}(x_{1}) + w_{1}f^{k−1}(x_{2}) + ⋯ + w_{k}f^{k−1}(x_{k}) = I_{k−1} } (∗)

where f^{m}(x) = x^{k}, and I_{m} = ∫_{a}^{b} f^{m}(x)dx.

• Romberg integration. Simpson’s rule turns out to be equivalent to a weighted sum of two applications of the Trapezoid Rule. Let I_{j} be the result of the trapezoid rule on 2^{j} + 1 points, and S_{j} the result of Simpson’s rule on the same grid. Then one can readily show that

S_{j+1} = \frac{4}{3}I_{j+1} − \frac{1}{3}I_{j}

and a careful error analysis shows that error terms in the trapezoid rules partially cancel, leaving an error that depends only on f^{(4)}(x)—in other words, even though Simpson’s involves only a quadratic approximation to f(x), it yields results that are exact up to cubics! Romberg integration methods are general methods for taking weighted sums of trapezoid rules to obtain cancellation of higher order terms like this.

2.2 Gaussian quadrature

\underline{Newton-Cotes} rules for integrating polynomials up to order k − 1 exactly are based on solving the linear equations (∗) for weights w_{1}, …, w_{k} at k equally spaced grid points x_{1}, …, x_{k},

∫_{a}^{b} f(x)dx = ∑_{i=0}^{n} w_{i}f(x_{i})

and then extending by adding the integral approximations over adjacent groups of k points.

\underline{Gaussian quadrature} methods allow not only the weights w_{1}, …, w_{k} to to be freely chosen, but also the locations x_{1}, …, x_{k} to be freely chosen. We can then obtain integration rules that are exact for polynomials up to order 2k − 1.

Gaussian quadrature formulae can also be obtained by solving the system (∗), but now (∗) is a nonlinear system. The nonlinear systems become ill conditioned quickly as k grows, so that either (a) very highprecision arithmetic must be used; or (b) more accurate methods employing stable recurrences relating othrogonal polynomials to one another are used.

Usually the problem is generalized to deal with integrals of the form

∫_{a}^{b} f(x)w(x)dx

for fixed weight function w(x) (which is equivalent to the density of a random variable with all moments finite).

Polynomials are chosen for the functions f^{j}(x) in (∗) that are orthogonal with respect to the w-inner product:

⟨f^{j_{1}}, f^{j_{2}}⟩_{w} = ∫_{a}^{b} f^{j_{1}}(x)f^{j_{2}}(x)w(x)dx = \{ \begin{matrix}0, j_{1} ≠ j_{2} \\ c^{j_{1}}, j_{1} = j_{2} \end{matrix}

Sometimes c^{j} ≡ 1 ∀j is desired (orthonormal polynomials) and sometimes c^{j} is allowed to be free and instead a standardization like f^{j}(1) ≡ 1 ∀j is chosen.

The solution of (∗) jointly for w_{i} and x_{i} is clearly tied up with the mathematics of the Hilbert space L_{w}^{2} of all functions f with finite integral ⟨f, f⟩_{w}. The main result is as follows:

Let w(x) be a density on (a, b) for a random variable X with all moments finite, and let f^{j}(x), j = 1, …, k be a set of orthonormal polynomials with respect to ⟨·, ·⟩_{w}, and let c_{j} be the leading coefficient of each f^{j}(x). Then:

• The positions a < x_{1} < ⋯ < x_{k} < b given by the roots of f^{k}(x), and the weights w_{j} > 0 given by

w_{j} = −\frac{c_{k+1}}{c_{k}} \frac{1}{f^{k+1}(x_{j})[f^{k}]′(x_{j})}

sastify

∫_{a}^{b} f(x)w(x)dx = ∑_{i=0}^{n} w_{i}f(x_{i})

whenever f(x) is a polynomial of degree 2k − 1 or less; and

• If f(x) is any 2k-times differentiable function,

∫_{a}^{b} f(x)w(x)dx − ∑_{i=0}^{n} w_{i}f(x_{i}) = \frac{f^{(2k)}(x^{∗})}{(2k)!c_{k}^{2}}, , , , x ∈ (a, b)

A very common case is when (a, b) = (−1, 1) and w(x) ∝ 1, and we take the standardization f^{k}(1) ≡ 1. In that case the orthogonal polynomials are the Legendre polynomials, and the quadrature method is known as Gauss-Legendre or just Gaussian quadrature.

Tables of x_{j}’s and w_{j}’s can be found in Abramowitz and Stegun (1964), but it is almost always better to use a “canned” routine rather than rederive these or program them by cutting and pasting from a table. It is almost never worthwhile, computationally or from an error-introducing point of view, to generate the x_{j}’s and w_{j}’s.

Automatic/extension forms of Gaussian quadrature rules are harder to come by than for Newton-Cotes methods, because the x_{1}, …, x_{k} do not form a subset of locations for some larger k, in general. The most commonly used approximate extensions are for Gauss-Legendre quadrature; they are due to Kronrod, and they are implemented in the Splus (and R) function integrate().

2.3 Numerical integration to compute posterior probabilities

Suppose we have obtained an interval [a, b] which we hope is “95% credible” for θ, i.e. P {θ ∈ [a, b]|y} ≈ 0.95. How can we verify this? The integrals in equation (1) are now of the form

p(θ ∈ [a, b]|y) = \frac{∫_{a}^{b} L(y|θ)π(θ)dθ}{∫_{−∞}^{∞} L(y|θ)π(θ)dθ}

These integrals present two immediate computational problems:

1. Computers don’t deal well with large numbers like ±∞;

2. As n grows, L(y|θ)π(θ) = ∏_{i=1}^{n} L_{i}(y_{i}|θ)π(θ) is likely to become quite small; computers don’t deal well with small numbers either.

For Problem 1. In a few problems, the support, or range of θ for which L(y|θ)π(θ) > 0, will actually be a finite interval, in which case 1. is not a problem. More likely the support will be (0, ∞), (1, ∞), or (−∞, ∞). If so, then we can try to use the change of variables formula

∫_{A}^{B} f(θ)dθ = ∫_{u(A)}^{u(B)} \frac{f(u^{−1}(t))}{u′(u^{−1}(t))}dt

where f(θ) ≡ L(y|θ)π(θ), and A or B may be infinite. We should choose the transformation u(θ) with two goals in mind:

• Make (u(A), u(B)) a nice interval. For example,

Transformation Some u(θ)’s
( 0,∞) → (0,1): e^{−θ}; θ/(1 + θ)
( 1,∞) → (0,1): 1/x
(−∞,∞) → (0,1): \frac{2}{π} tan^{−1}(θ); ; e^{θ}/(1 + e^{θ}); ; other CDF’s.

• Make g(t) = \frac{f(u^{−1}(t))}{u′(u^{−1}(t))} nearly equal to some easy-to-integrate function (constant, polynomial, etc.)

This is difficult; for example to ask for g(t) = constant is basically to ask for f/u′ = constant, i.e.

u(θ) = ∫ f′, dθ

The next best thing to do is go for computational convenience, or match tails, center, or spread of f with u′.

An alternative to the transformation method is just to graph (something proportional to) f(θ) choose A′ and B′ from the graph, so that the integral from A′ to B′ is a good approximation to the integral from −∞ to ∞.

For Problem 2. If L(y|θ)π(θ) is small, the basic strategy for fixing this is to work with c · L(y|θ)π(θ) instead. The idea is that c will cancel out of the numerator and denominator integrals in equation (1). Choosing c can be difficult. Some first steps include:

• Graph c · L(y|θ)π(θ) as a function of θ for various trial values of c and see what happens

• Divide by L(y|θ̂)π(θ̂) where θ̂ is the MLE or posterior mode (probably you want to work with the difference of the logs and then exponentiate).

3 Sampling to produce pictures of posterior densities; Monte Carlo integration

This section gives the background for sketches two simple sampling methods to produce samples from the posterior distribution in a Bayesian inference problem. Some books to consider for further details about Monte Carlo simulation are:

• Devroye, L. (1986). \underline{Nonuniform random variate generation.} New York: Springer-Verlag.

• Hammersley, I. M. and Handscomb, D. C. (1964). \underline{Monte Carlo methods.} London: Metheun.

• Ripley, B. D. (1987). \underline{Stochastic simulation.} New York: Wiley.

• Gelman, A., Carlin, J. Stern, H. and Rubin, D. (1996?) Bayesian data analysis. New York: Chapman-Hall.

The connections with Bayesian posterior sampling come from the 1992 article by Smith and Gelfand that I am also handing out this week.

3.1 Sampling Generalities

Monte Carlo sampling typically has either (or both) of two goals:

1. To generate a sample from a distribution so that it may be plotted as a histogram, density estimate, etc.

2. To evaluate an integral, such as the mean, variance, etc., of a distribution.

For either of these goals the approach is the same:

Produce a sample that is (approximately) an i.i.d. sample from the distribution of interest.

Most programming languages used in statistical calculations (including, of course, SPLUS) have built-in pseudo-random number generators that generate observations that appear to be i.i.d. Unif[0, 1]. In this section we’ll list some common ways to generate observations from other distributions using these.

As a matter of notation, let

f(x) = the density of the target distribution
F (x) = ∫_{−∞}^{x} f(t)dt = P [X ≤ x] = the CDF for f(x)
g(y) = the density of a convenient distribution
G(y) = ∫_{−∞}^{y} g(t)dt = P [Y ≤ x] = the CDF for g(y)

Then a generic way to satisfy each of the goals above is:

1. Generate x_{1}, …, x_{n} as i.i.d. observations from f(x).

2. With the data generated in 1., compute

\frac{1}{n} ∑_{i=1}^{n} h(x_{i}) ≈ ∫_{−∞}^{∞} h(x)f(x)dx ≡ E_{f}[h(X)]

Probability Integral Transform You can verify with calculus the following two identities.

• If X is distributed following f(x), then U = F (X) is Unif[0, 1].

• If U is Unif[0, 1], then F^{−1}(U) is distributed following f(x).

If F (x) can be obtained and inverted in closed form, this provides a way to generate from X. Relatively few distributions are easy to sample from in this way, however.

3.1.1 Clever Transformations

One often learns about these in math stat courses as “applications” of the “change of variables” formula in integral calculus. They actually turn out to be useful and convenient in some cases, though most of the distributions for which there are easy transformations are already coded into standard random number function (e.g. in SPLUS).

The basic idea is that if h(x) is a strictly increasing function, then

Y ∼ g(y) ⇒ X = h^{−1}(Y) ∼ g(h(x))h′(x)

This is easily seen from the calculation

\frac{d}{dx}P [X ≤ x] = \frac{d}{dx}P [h^{−1}(Y) ≤ x]
= \frac{d}{dx}P [Y ≤ h(x)] = \frac{d}{dx}G(h(x))
= g(h(x))h′(x)

(there are many details, such as what to do if h(x) is strictly decreasing or even noninvertible but this is the basic idea.)

If you want to sample from f(x) then the goal is to find h(x) so that f(x) = g(h(x))h′(x). Sometimes one finds h by accident, sometimes by experience, and sometimes by explicitly solving this equation. There is a similar but more complicated formula when x is multivariate. Here are some typical examples of transformations that lead to familiar distributions, in both the univariate and multivariate cases:

• Y ∼ Unif[0, 1], h^{−1}(y) = − log y ⇒ X = h^{−1}(Y) ∼ Exp(1).

• Y ∼ Unif[0, 1] ⇒ X = tan[π(Y − 1/2)] ∼ Cauchy (this is also an instance of the probability integral transform).

• Y_{1}, Y_{2} ∼ independent Unif[0, 1]’s ⇒ X_{1} = cos(2πY_{1})\sqrt{−2 log(Y_{2})} and X_{2} = sin(2πY_{1})\sqrt{−2 log(Y_{2})} are independent N(0, 1)’s.

• Y ∼ N(0, 1) ⇒ Y^{2} ∼ Chi-squared, df = 1.

• Y_{1}, …, Y_{n} ∼ independent N(0, 1)’s ⇒ X = Y_{1}^{2} + Y_{2}^{2} + ⋯ + Y_{n}^{2} ∼ Chi-squared, df = n.

3.1.2 Rejection Method

More frequently the cdf F (x) will not be readily computable nor invertable, and there will be no obvious transformation to produce X ∼ F (x) from uniform or other simpler random variables. The rejection method is a useful, generic approach to sampling from such problematic F (x)’s.

The Rejection Method. Suppose f(x) is \underline{proportional to} the density we \underline{want} to sample from; and g(x) is a density we \underline{can} sample from. Also, suppose

f(x) ≤ a · g(x) for all x

Then the following procedure produces observations from f(x):

• Generate Y from g(x);

• Generate U from Unif[0, a · g(x)];

• If U < f(Y), accept Y as an observation X from f(x); otherwise reject Y and start over.

Useful exercises:

• Use calculus in 2D, to show that P [Y ≤ x|Y is accepted] = F (x), justifying the method.

• Show that you can generate X ∼ N(0, 1) by applying the rejection method to Cauchy observations Y = tan(π(U − 1/2)). What is (the \underline{best}) a?

• Show (analytically) that it does \underline{not} work to generate X ∼ Cauchy by applying the rejection method to N(0, 1) observations Y. What goes wrong???

The fact that one does not need to know the “normalizing constant” for f(x) is a big advantage, especially in Bayesian problems, as we shall see below.

3.1.3 Control Variates

The basic idea of control variates is that if you want to estimate an integral with an i.i.d. sum,

E_{f}[h(X)] = ∫ h(x)f(x)dx ≈ Ê_{n} ≡ \frac{1}{n} ∑_{i=1}^{n} f(X_{i})

and you can arrange to get another i.i.d sum that estimates 0,

0 ≈ Ẑ_{n}

\underline{and} Cov(Ê_{n}, Ẑ_{n}) < 0, then Ê_{n} + Ẑ_{n} ≈ E_{f}[h(X)] still works, and

Var(Ê_{n} + Ẑ_{n}) = Var(Ê_{n}) + 2Cov(Ê_{n}, Ẑ_{n}) + Var(Ẑ_{n}) < Var(Ê_{n})

as long as 2|Cov(Ê_{n}, Ẑ_{n})| > Var(Ẑ_{n}).

A similar trick works with ratios Ê_{n}/Ẑ_{n} where now Ẑ_{n} ≈ 1 and Cov(Ê_{n}, Ẑ_{n}) > 0; an example of this “ratio method” will be given in the next subsection.

3.1.4 Importance Sampling and Sampling-Importance-Resampling (SIR)

Importance sampling Importance sampling is natural when one has to calculate E_{f}[h(X)] where X ∼ f(x) but it is easier to sample from g(x). There are two versions of importance sampling.

\underline{Simple importance sampling.}Since

E_{f}[h(X)] = ∫ h(x)f(x)dx = ∫ \frac{h(x)f(x)}{g(x)} · g(x)dx = E_{g}[h(Y)f(Y)/g(Y)]

we can estimate E_{f}[h(X)] as follows:

• Generate Y_{1}, …, Y_{n} i.i.d. from g(y);

• Estimate E_{f}[h(X)] ≈ \frac{1}{n} ∑_{1}^{n} h(Y_{i})f(Y_{i})/g(Y_{i}).

\underline{Fancy importance sampling.} Sometimes, especially in Bayesian problems, it is useful to be able to do importance sampling without computing the normalizing constant for f(x). We noted above that

Ĩ_{n} = \frac{1}{n} ∑_{1}^{n} h(Y_{i})f(Y_{i})/g(Y_{i}) → E_{f}[h(X)]

It similarly follows that

Õ_{n} = \frac{1}{n} ∑_{1}^{n} 1 · f(Y_{i})/g(Y_{i}) → 1

so a fancier version of importance sampling that is sometimes used is

Î_{n} ≡ \frac{Ĩ_{n}}{Õ_{n}} = \frac{∑_{1}^{n} h(Y_{i})f(Y_{i})/g(Y_{i})}{∑_{1}^{n} f(Y_{i})/g(Y_{i})} → E_{f}[h(X)]

This version has two possible advantages.

• Since f(x) appears in both the numerator and denominator, any normalizing constant would cancel, so it is only necessary to know f(x) up to a(ny) multiplicative constant. Similarly, g(x) need only be known up to a multiplicative constant.

• We can see that the variance of the simpler version is

Var(Ĩ_{n}) = \frac{1}{n}Var_{g}[h(Y)f(Y)/g(Y)]

and, applying the multivariate delta method to the function g(x, y) = x/y, where x = \frac{1}{n} Ĩ_{n} and y = \frac{1}{n} Õ_{n}, the variance of the fancier version is approximately

Var(Î_{n}) = Var(Ĩ_{n}) − 2E_{f}[h(X)]Cov(Ĩ_{n}, Õ_{n}) + (E_{f}[h(X)])^{2}Var(Õ_{n})

By making Cov(Ĩ_{n}, Õ_{n}) large and Var(Õ_{n}) small, we can possibly make the MC SE of the fancier version of importance sampling less than the MC SE of the simpler version.

Sampling Importance Resampling (SIR) This is a method of using importance sampling to generate a sample from a new distribution, not just estimate an integral. It also goes by the name “weighted bootstrap”.

Note that, in the fancier version of importance sampling, if we take h(x) to be

h(x) = 1_{{x≤x_{0}}}

then, since h(Y_{i}) = 0 if Y_{i} > x_{0} and = 1 if Y_{i} ≤ x_{0},

Î_{n} = \frac{∑_{1}^{n} h(Y_{i})f(Y_{i})/g(Y_{i})}{∑_{1}^{n} f(Y_{i})/g(Y_{i})} = ∑_{i: Y_{i}≤x_{0}} q_{i}

where q_{i} = [f(Y_{i})/g(Y_{i})]/ ∑_{1}^{n} f(Y_{j})/g(Y_{j}). But also,

Î_{n} = ∑_{i: Y_{i}≤x_{0}} q_{i} ≈ E_{f}[1_{{X≤x_{0}}}] = P [X ≤ x_{0}] = F (x_{0}).

Now consider the discrete random variable T that takes on the value Y_{i} with probability q_{i}; the cdf of T is

P [T ≤ x_{0}] = ∑_{i: Y_{i}≤x_{0}} q_{i} = Î_{n} ≈ P [X ≤ x_{0}]

This leads to the following method of modifying a random sample Y_{1}, …, Y_{n} from g(y) so that it looks like a random sample X_{1}, …, X_{n} from f(x).

SIR. Suppose f(x) is \underline{proportional to} the density we \underline{want} to sample from; and g(x) is a density we \underline{can} sample from. Then the following procedure produces observations that are \underline{approximately} distributed according to f(x):

• Generate Y_{1}, …, Y_{M} i.i.d. from g(x);

• Compute q_{1}, …, q_{M} as q_{i} = [f(Y_{i})/g(Y_{i})]/ ∑_{1}^{n} f(Y_{j})/g(Y_{j});

• Sample X_{1} …, X_{n} i.i.d. from the Y ’s with probabilities the q’s.

Once again, we really only need to know f(x) up to a multiplicative constant.

A cautionary note: Usually you need M ≫ n (like, at least M > 10n) to make the distribution of the X’s reasonably close to f(x).

3.2 Exploring posterior distributions in Bayesian inference, using Monte Carlo

Given a (possibly multidimensional) parameter θ with prior distribution π(θ) and likelihood L(y|θ) = ∏_{i} L_{i}(y_{i}|θ), Bayesian inference is based on functionals (i.e. numerical functions) of the posterior distribution

p(θ|y) = \frac{L(y|θ)π(θ)}{∫ L(y|θ)π(θ)dθ}

Common functionals are

• The posterior probability of an interval [a, b], ∫ 1_{{θ∈[a,b]}}p(θ|y)dθ;

• The posterior mean E[Θ|y] = ∫ θp(θ|y)dθ;

• The posterior variance Var(Θ|y) = ∫ (θ − E[Θ|y])^{2}p(θ|y)dθ;

• The “entire” posterior density p(θ|y).

A direct attack on the first three quantities requires at the very least two integrals:

E[h(θ)|y] = ∫ h(θ)p(θ|y)dθ = \frac{∫ h(θ)L(y|θ)π(θ)dθ}{∫ L(y|θ)π(θ)dθ}  (2)

In this section we show how rejection and importance sampling provide basic tools for estimtating the integral (1) and/or obtaining a sample from the full posterior distribution p(θ|y).

3.2.1 Rejection Method

To link up with the general rejection method, let

f(θ) = L(y|θ)π(θ), and
g(θ) = π(θ).

In many Bayesian problem s g(θ) will be relatively simple to sample from, e.g. it may be a product of univariate densities like normals or gammas or uniforms, that SPLUS has built-in functions to generate samples from, so that each coordinate of θ may be sampled independently from its own SPLUS generator. p(θ|y) is usually hard to sample from, and we typically only know f(θ) ∝ p(θ|y).

For the rejection method, we need

f(θ) ≤ a · g(θ), or
L(y|θ)π(θ) ≤ a · π(θ)

This is easily accomplished by taking a = L(y|θ̂) where θ̂ is the MLE for the problem. Thus the rejection method becomes:

Rejection for Bayes.

• Generate θ from π(θ);

• Generate U from Unif[0, L(y|θ̂)π(θ)];
• If U < L(y|θ)π(θ), accept θ as an observation from p(θ|y);

otherwise reject θ and start over.

Note that π(θ) can be dropped from steps 2 and 3, if we generate U ∼ Unif[0, 1] and check U < L(y|θ)/L(y|θ̂). This is good.

Also note that if the prior π(θ) disagrees radically with the likelihood L(y|θ) about where θ lives (compare the modes, for example), then this is going to be a very slow way to generate posterior observations. This is bad.

3.2.2 Importance Sampling and SIR

A useful question to apply Importance Sampling and SIR to in Bayesian inference is, “What if I had chosen a different prior distribution?”

Assume that θ_{1}, …, θ_{M} are a random sample from p(θ|y) ∝ L(y|θ)π(θ), and we want switch the prior from π(θ) to π^{∗}(θ).

Now if we let

f(θ) = L(y|θ)π^{∗}(θ), and
g(θ) = L(y|θ)π(θ)

we see that the “importance weights” q_{i} become

q_{i} = [f(θ_{i})/g(θ_{i})]/ ∑_{1}^{n} f(θ_{j})/g(θ_{j}) = [π^{∗}(θ_{i})/π(θ_{i})]/ ∑_{1}^{n} π^{∗}(θ_{j})/π(θ_{j})

Thus, posterior expected values for the new prior π^{∗} can be approximated using importance sampling as

E_{π^{∗}}[h(θ)|y] ≈ ∑_{1}^{M} h(θ_{i}) q_{i}, for i.i.d. θ_{i} ∼ p(θ|y) ∝ L(y|θ)π(θ)

and an SIR sample (approximately) from the new posterior can be obtained by

Sample θ_{1}^{∗}, …, θ_{n}^{∗} as i.i.d. from θ_{1}, …, θ_{M}, with probabilities q_{1}, … q_{M}.