Sampling methods Hans-Peter Helfrich University of Bonn - - PowerPoint PPT Presentation

Sampling methods

Hans-Peter Helfrich

University of Bonn Theodor-Brinkmann-Graduate School

June 5, 2013

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 1 / 19

Overview

1

Introduction

2

Acceptance Rejection sampling

3

Gibbs’ sampling

4

Markov Chain Monte Carlo algorithm

5

References

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 2 / 19

Introduction

Bayes’ theorem

Statistical inference can be done via Bayes’ theorem. The posterior distribution is given by p(휃∣y) ∝ L(y∣휃)p(휃) The posterior density function contains all statistical information for providing mean values, medians, and measures of uncertainty.

Sampling methods

In general, the normalizing factor cannot be explicitly calculated. That can be done by sampling methods that give samples with the posterior density

• distribution. We consider

Acceptance Rejection sampling Gibbs’ sampling Markov Chain Monte Carlo Method

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 3 / 19

Sampling

Methods

There are many methods for sampling. We mainly consider Acceptance Rejection sampling Gibbs’ sampling is applied in WinBUGS [Lunn et al., 2000] and the

• pen source counterpart OpenBUGS, see The Bugs Project. For

small to medium sized problems, both programs are very convenient for parameter estimation. The methods can be called from R, Matlab and other statistical software. For large problems with many parameters and complex physical models, it may be advantageous to build your own Markov Chain Monte Carlo implementation, see, e.g. [Van Oijen et al., 2005]. Here prediction in forest models is done by a complex model, which utilizes this method.

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 4 / 19

Acceptance Rejection sampling

Algorithm

We want to sample from a density distribution f (x). We choose a distribution g(x), where g need not to be a proper distribution. We denote M the maximum of g(x). We may also choose the non-informative g(x) = 1. Sample x from g(x) and u from the uniform distribution over the unit interval. Check if u ≤ f (x) Mg(x)

If this holds, accept x as a realization of f (x) Otherwise, reject the value of x and repeat the sampling step

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 5 / 19

Gibbs’ sampling

Assumptions

We want to draw samples from a n-dimensional distribution. We consider the case n = 2 and denote the density function by f (x, y).

The Gibbs’ sampler

This method is also called alternating conditional sampling. At each iteration k the Gibbs’ sampler cycles through both components xk+1 is sampled from f (x∣yk) yk+1 is sampled from f (y∣xk+1)

Conditional distributions

The conditional distributions are given by f (x∣y) ∝ f (x, y), f (y∣x) ∝ f (x, y)

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 6 / 19

Sampling a two-dimensional density distribution

We want to sample from a two-dimensional distribution with the density f (x, y) ∝ exp ( −2 3(x2 + 4y2 − 2xy) )

Conditional densities

f (x∣y) ∝ exp ( −2 3(x2 + 4y2 − 2xy) ) ∝ exp ( −2 3(x − y)2 ) f (y∣x) ∝ exp ( −2 3(x2 + 4y2 − 2xy) ) ∝ exp ( −2 3(2y − x 2)2 )

Conditional densities

f (x∣y) = N(x∣y, 3 4), f (y∣x) = N(y∣x 4, 3 16)

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 7 / 19

Example

Gibbs’ sampler

Start with arbitrary x0, y0, then iterate xk+1 is drawn from 풩(yk, 3 4) yk+1 is drawn from 풩(xk+1 4 , 3 16) In order to suppress the dependence on the starting values, the first, say 1000 or 10000 iterations should be discarded.

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 8 / 19

Simulation in R

N <- 2000 x <- rep(0,N) y <- rep(0,N) for(k in (1:(N-1))) { x[k+1] <- rnorm(1, y[k], sqrt(3/4)) y[k+1] <- rnorm(1, x[k + 1]/4, sqrt(3/16)) } print(sd(x)) print(sd(y)) print(cor(x,y))

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 9 / 19

Figure

We get a picture of the samples as pdf-file by the commands pdf("gibbs.pdf") plot(x,y, xlim = c(-3,3), ylim = c(-3,3), col ="red") dev.off()

• ●
• ●
• ●
• ●
• ●
• ●
• −3

−2 −1 1 2 3 −3 −2 −1 1 2 3 x y

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 10 / 19

General idea of MCMC

Generating samples by jumping from location x to y

휋(x), 휋(y) probability densities at locations x, y Jumping from x to y with probability density q(x, y) = q(y, x) x → y accept y with probability 휋(y)/휋(x) y → x accept x x → y 휋(x) 휋(y)

Goal

Number of samples at y Number of samples at x = 휋(y) 휋(x)

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 11 / 19

Model

Nonlinear compartment model

We want to estimate the parameters a, k1, k2 in our nonlinear compartment model c(t) = ak1 (k1 − k2)(exp(−k2t) − exp(−k1t)) (1) for given measurements c(t1), c(t2), . . . at times t1, t2, . . . .

Statistical model

In order to keep k1 > k2, we introduce d = k1 − k2 > 0 as parameter. We have measurements of the concentrations c(tk) at times t1, t2, . . . , tN. cj = c(tj, a, k1, d) + 휖j, 휖j ∈ 풩(0, 휎2)

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 12 / 19

Statistical model in OpenBUGS

Model

For the parameters a, k1, d and the precision 휏 of the the error term 휖i, we choose non-informative prior density distributions. regr <- function() { k2 <- k1 - d for (j in 1:N) { conc[j] ˜ dnorm (theta[j], tau) theta[j] <- a*k1/(k1-k2)*(exp(-k2*x[j]) - exp(-k1*x[j])) } a ˜ dunif(5,100) k1 ˜ dunif(0.3,10) d ˜ dunif(0,2) tau ˜ dunif(0,10) # non-informative prior on tau }

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 13 / 19

Bayesian modelling with BUGS

The BUGS code given before may be used in several environments: WinBUGS may be used as stand-alone program. You may get it free of charge from http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml. The RBugs package [Gelman et al., 2004] http://www.stat.columbia.edu/ gelman/bugsR/software.pdf may be used for calling WinBUGS within R. All necessary steps can be done within the graphical user interface RGUI of R. The R2WinBUGS package [Sturtz et al., 2005] provides convenient functions to call WinBUGS or the open source counterpart OpenBUGs from R MATBUGS http://code.google.com/p/matbugs/ gives a Matlab interface for WinBUGS.

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 14 / 19

Parameter estimates

Estimates by nonlinear least squares

Estimate Std. Error t value Pr(>|t|) a 48.821 29.822 1.64 0.16 k1 1.540 1.292 1.19 0.29 d 0.935 1.789 0.52 0.62 Residual standard error: 3.04 on 5 degrees of freedom

Estimates by BUGS

mean sd 2.5% 25% 50% 75% 97.5% a 50.33 11.84 31.45 40.84 48.97 59.30 73.12 k1 1.71 0.65 0.93 1.22 1.54 2.01 3.45 d 1.09 0.82 0.03 0.45 0.94 1.53 3.15 tau 0.16 0.09 0.04 0.10 0.15 0.21 0.37 Uncertainty is visualized by the density distributions

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 15 / 19

Output of OpenBUGS

80% interval for each chain R−hat 50 50 100 100 1 1.5 2+ 1 1.5 2+ 1 1.5 2+ 1 1.5 2+ 1 1.5 2+ 1 1.5 2+ a

• k1
• d
• tau
• medians and 80% intervals

a

20 40 60 80

• k1

1 2 3

• d

1 2 3

• tau

0.1 0.2 0.3

• deviance

35 40 45

• Bugs model at "theanin.txt", fit using OpenBUGS, 4 chains, each with 20000 iterations (first 10000 discarded)

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 16 / 19

Distribution of the parameter estimates

Parameter a

adist Density 20 40 60 80 100 0.000 0.015 0.030

Parameter k1

k1dist Density 1 2 3 4 0.0 0.4 0.8

Parameter d

ddist Density 1 2 3 4 0.0 0.2 0.4 0.6

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 17 / 19

Fitting of the data

• 5

10 15 20 5 10 15 20 25 30

Fitting a compartment model

Time Concentration Openbugs Nonlinear least squares H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 18 / 19

References

Gelman, A. B., Carlin, J. S., Stern, H. S., and Rubin, D. B. (2004). Example of Computation in R and Bugs, chapter Appendix C, pages 591 – 609. Chapman & Hall. Lunn, D. J., Thomas, A., Best, N., and Spiegelhalter, D. (2000). Winbugs - a Bayesian modelling framework: Concepts, structure, and extensibility. Statistics and Computing (2000), 10:325–337. Sturtz, S., , Ligges, U., and Gelman, A. (2005). R2WinBUGS: A package for running WinBUGS from R. Journal of Statistical Software, 12:1–16. Van Oijen, M., Rougier, J., and Smith, R. (2005). Bayesian calibration of process-based forest models: bridging the gap between models and data. Tree Physiology, 25(7):915–927.

H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 19 / 19