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 Introduction 1 Acceptance Rejection sampling 2 Gibbs’ sampling 3 Markov Chain Monte Carlo algorithm 4 References 5 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 open 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 f ( x ) u ≤ 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 f ( x ∣ y k ) x k +1 is sampled from y k +1 is sampled from f ( y ∣ x k +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 ( − 2 ) 3( x 2 + 4 y 2 − 2 xy ) f ( x , y ) ∝ exp Conditional densities ( − 2 ) ( − 2 ) 3( x 2 + 4 y 2 − 2 xy ) 3( x − y ) 2 f ( x ∣ y ) ∝ exp ∝ exp ( − 2 ) ( − 2 3(2 y − x ) 3( x 2 + 4 y 2 − 2 xy ) 2) 2 f ( y ∣ x ) ∝ ∝ exp exp Conditional densities f ( x ∣ y ) = N( x ∣ y , 3 f ( y ∣ x ) = N( y ∣ x 4 , 3 4) , 16) H.-P. Helfrich (University of BonnTheodor-Brinkmann-Graduate School) Sampling methods June 5, 2013 7 / 19

Example Gibbs’ sampler Start with arbitrary x 0 , y 0 , then iterate 풩 ( y k , 3 x k +1 is drawn from 4) 풩 ( x k +1 , 3 y k +1 is drawn from 16) 4 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