Color test I just want to be sure that all my favourite colors are being displayed correctly on this new device. If not I’ll modify them. R´ emi Bardenet Bayesian numerical methods 1 / 51 How to unefficiently solve any problem A tutorial on Bayesian numerical methods R´emi Bardenet LAL, LRI, Univ. Paris-Sud XI May 30th, 2012 R´ emi Bardenet Bayesian numerical methods 2 / 51 Summary 1 A very generic problem 2 First sampling methods 3 MCMC algorithms 4 A taste of a monster MCMC sampler for Auger R´ emi Bardenet Bayesian numerical methods 3 / 51 Summary 1 A very generic problem 2 First sampling methods 3 MCMC algorithms 4 A taste of a monster MCMC sampler for Auger R´ emi Bardenet Bayesian numerical methods A very generic problem 4 / 51 Today’s talk: sampling for Bayesian Inference When I a model p(data|θ) of an experiment has been written, I a prior p(θ) has been set on the parameters, then by Bayes’ theorem: π(θ) := p(θ|data) = R R´ emi Bardenet Bayesian numerical methods p(data|θ)p(θ) . Θ p(data|θ)p(θ)dθ A very generic problem 5 / 51 Today’s talk: sampling for Bayesian Inference When I a model p(data|θ) of an experiment has been written, I a prior p(θ) has been set on the parameters, then by Bayes’ theorem: π(θ) := p(θ|data) = R R´ emi Bardenet Bayesian numerical methods p(data|θ)p(θ) . Θ p(data|θ)p(θ)dθ A very generic problem 5 / 51 Today’s talk: sampling for Bayesian Inference When I a model p(data|θ) of an experiment has been written, I a prior p(θ) has been set on the parameters, then by Bayes’ theorem: π(θ) := p(θ|data) = R R´ emi Bardenet Bayesian numerical methods p(data|θ)p(θ) . Θ p(data|θ)p(θ)dθ A very generic problem 5 / 51 The generative process of an air shower R´ emi Bardenet Bayesian numerical methods A very generic problem 6 / 51 0.2 Looking at Auger tank signals 0.2 0.0 0 50 100 150 tbin [25 ns] 200 250 300 0.1 0.1 0.2 0.3 !V 0.2 0.3 !V 0.2 0.3 !V Event 1524469 - E = 10.6 [EeV] / ! = 39.6° - Station 612 at 1684 m 2.2 2.0 2 1.8 1.4 1 dN d! V VEM-peak 1.6 1.2 1.0 0.8 0.6 0.2 0.4 0.2 0.0 0 50 100 150 tbin [25 ns] 200 250 300 0.1 0.1 Event 1524469 - E = 10.6 [EeV] / ! = 39.6° - Station 622 at 2031 m 1.8 1.6 2 1.2 1 dN d! V VEM-peak 1.4 1.0 0.8 0.6 0.4 0.2 0.2 0.0 R´ emi Bardenet 0 50 100 150 tbin [25 ns] 200 250 300 0.1 0.1 F IG . 5.17 - Traces FADC sautsproblem pour les cinq stations Bayesian numerical methodset distributionsA des very generic 7 / 51 Modelling tank signals: the single muon case Figure: Muonic time response model pτ,td 0.015 p 0.010 0.005 0.000 -20 0 20 40 t @nsD 60 80 100 Mean number of Photo-electrons per bin Z ti n¯i (Aµ , tµ ) = Aµ pτ,td (t − tµ )dt. ti−1 R´ emi Bardenet Bayesian numerical methods A very generic problem 8 / 51 50 45 40 35 ðPE 30 25 20 15 10 5 0 0 100 200 300 400 500 600 t @nsD ni Poisson with mean n¯i (Aµ , tµ ) = Nµ X n¯i (Aµ j , tµ j ), j=1 R´ emi Bardenet Bayesian numerical methods A very generic problem 9 / 51 So what is what in the cas Nµ = 1 ? I Parameters to infer are θ = (tµ , Aµ ). I The likelihood is fixed by the model: p(data|θ) = N bins Y i=1 I Choose e.g. a uniform prior: p(θ) = I n¯i (θ)ni −¯ni (θ) e . ni ! 1 1 χ (tµ ) χ[0,D] (Aµ ), C D [0,C ] Then the posterior reads p(θ|data) ∝ χ[0,C ] (tµ ) χ[0,D] (Aµ ) N bins Y n¯i (θ)ni e −¯ni (θ) . i=1 R´ emi Bardenet Bayesian numerical methods A very generic problem 10 / 51 Bayesian inference requires computing integrals I MEP estimate (aka Bayes’): compute Z ˆ θMEP := θp(θ|data)dθ, Θ I credible intervals: find I such that Z p(θ|data)dθ ≥ 1 − α. I I Bayes’ factors: require evidence computations Z p(data|M) = p(data|θ, M)p(θ|M)dθ. The MAP estimate θMAP := arg max p(θ|data) with “Hessian” credible interval does not require integrals. R´ emi Bardenet Bayesian numerical methods A very generic problem 11 / 51 Averaging whenever necessary Most Bayesian inference tasks require integrals wrt. the posterior Z h(θ)p(θ|data)dθ. The Monte Carlo principle Z N 1 X I := h(x)π(x)dx ≈ h(xi ) =: ˆIN , N i=1 where X1:N ∼ π I I i.i.d. ˆIN is unbiased: EˆIN = I . √ Error bars shrink at speed 1/ N: Var(X ) . Var(ˆIN ) = N R´ emi Bardenet Bayesian numerical methods A very generic problem 12 / 51 Averaging whenever necessary Most Bayesian inference tasks require integrals wrt. the posterior Z h(θ)p(θ|data)dθ. The Monte Carlo principle Z N 1 X I := h(x)π(x)dx ≈ h(xi ) =: ˆIN , N i=1 where X1:N ∼ π I I i.i.d. ˆIN is unbiased: EˆIN = I . √ Error bars shrink at speed 1/ N: Var(X ) . Var(ˆIN ) = N R´ emi Bardenet Bayesian numerical methods A very generic problem 12 / 51 Averaging whenever necessary Most Bayesian inference tasks require integrals wrt. the posterior Z h(θ)p(θ|data)dθ. The Monte Carlo principle Z N 1 X I := h(x)π(x)dx ≈ h(xi ) =: ˆIN , N i=1 where X1:N ∼ π I I i.i.d. ˆIN is unbiased: EˆIN = I . √ Error bars shrink at speed 1/ N: Var(X ) . Var(ˆIN ) = N R´ emi Bardenet Bayesian numerical methods A very generic problem 12 / 51 Pros & Cons Pros I Numerical integration (grid methods) is too costly and imprecise when d ≥ 6. I MC should concentrate the effort on places where π is big. I Many, many applications ! Cons One must be able to sample from π ! → Need for generic clever sampling methods ! R´ emi Bardenet Bayesian numerical methods A very generic problem 13 / 51 Pros & Cons Pros I Numerical integration (grid methods) is too costly and imprecise when d ≥ 6. I MC should concentrate the effort on places where π is big. I Many, many applications ! Cons One must be able to sample from π ! → Need for generic clever sampling methods ! R´ emi Bardenet Bayesian numerical methods A very generic problem 13 / 51 Pros & Cons Pros I Numerical integration (grid methods) is too costly and imprecise when d ≥ 6. I MC should concentrate the effort on places where π is big. I Many, many applications ! Cons One must be able to sample from π ! → Need for generic clever sampling methods ! R´ emi Bardenet Bayesian numerical methods A very generic problem 13 / 51 Pros & Cons Pros I Numerical integration (grid methods) is too costly and imprecise when d ≥ 6. I MC should concentrate the effort on places where π is big. I Many, many applications ! Cons One must be able to sample from π ! → Need for generic clever sampling methods ! R´ emi Bardenet Bayesian numerical methods A very generic problem 13 / 51 Non-Bayesian applications of sampling algorithms: Auger “A complete Monte Carlo hybrid simulation has been performed to study the trigger efficiency and the detector performance. The simulation sample consists of about 6000 proton and 3000 iron CORSIKA [19] showers” from “The exposure of the hybrid detector of the P. Auger Observatory”, the Auger Collab., Astroparticle Phys., 2010. R´ emi Bardenet Bayesian numerical methods A very generic problem 14 / 51 Non-Bayesian applications of sampling algorithms: EPOS “The difficulty with Monte Carlo generation of interaction configurations arises from the fact that the configuration space is huge and rather nontrivial [...] the only way to proceed amounts to employing dynamical MC methods” from “Parton-based Gribov-Regge theory”, Drescher et al., Phys.Rept., 2008. R´ emi Bardenet Bayesian numerical methods A very generic problem 15 / 51 Summary 1 A very generic problem 2 First sampling methods 3 MCMC algorithms 4 A taste of a monster MCMC sampler for Auger R´ emi Bardenet Bayesian numerical methods First sampling methods 16 / 51 The inverse function method Principle If U ∼ U(0,1) and F is a cdf, then F −1 (U) ∼ F . I It assumes we know how to sample from U(0, 1) ! I It needs a known and convenient cdf. R´ emi Bardenet Bayesian numerical methods First sampling methods 17 / 51 The inverse function method Principle If U ∼ U(0,1) and F is a cdf, then F −1 (U) ∼ F . I It assumes we know how to sample from U(0, 1) ! I It needs a known and convenient cdf. R´ emi Bardenet Bayesian numerical methods First sampling methods 17 / 51 Rejection Sampling 1/2 Assumptions I Target π is known up to a multiplicative constant, I an easy-to-sample distribution q s.t. π ≤ Mq is known. RejectionSampling π, q, M, N 1 for i ← 1 to N, 2 Sample θ(i) ∼ q and u ∼ U(0,1) . 3 Form the acceptance ratio ρ = 4 if u < ρ, then accept θ(i) . 5 else reject. R´ emi Bardenet Bayesian numerical methods π(θ(i) ) . Mq(θ(i) ) First sampling methods 18 / 51 Rejection Sampling 2/2 Remarks & Drawbacks I One needs to know good q and M, I Lots of wasted samples. R´ emi Bardenet Bayesian numerical methods First sampling methods 19 / 51 Importance Sampling 1/2 Assumptions & principle I Target π is fully known, I an easy-to-sample distribution q is known s.t. Supp(π) ⊂ Supp(q). Then Z N 1 X π(θ(i) ) h(θ(i) ) IˆN := → h(θ)π(θ)dθ, N q(θ(i) ) Θ θ(i) ∼ q i.i.d. i=1 ImportanceSampling π, q, N 1 Sample independent θ(i) ∼ q, i = 1..N, 2 Form the weights wi = 3 π is approximated by R´ emi Bardenet Bayesian numerical methods 1 N π(θ(i) ) . q(θ(i) ) PN i=1 wi δθ(i) . First sampling methods 20 / 51 Importance Sampling 2/2 I If |h| ≤ M, then M2 Var(IˆN ) ≤ N I Z (π − q)2 dθ + 1 , q thus q has to be close to π and have heavier tails than π. Further: the q achieving the smallest variance is q=R |h|π . |h|πdθ I Better and achievable: For smaller variance, or if π is known up to a constant, use PN h(θ(i) )π(θ(i) )/q(θ(i) ) ˜ IN := i=1 (biased). PN (i) (i) i=1 π(θ )/q(θ ) I One needs to know a good, heavy-tailed q. I Adaptive strategies for tuning q are possible [WKB+ 09]. R´ emi Bardenet Bayesian numerical methods First sampling methods 21 / 51 Importance Sampling 2/2 I If |h| ≤ M, then M2 Var(IˆN ) ≤ N I Z (π − q)2 dθ + 1 , q thus q has to be close to π and have heavier tails than π. Further: the q achieving the smallest variance is q=R |h|π . |h|πdθ I Better and achievable: For smaller variance, or if π is known up to a constant, use PN h(θ(i) )π(θ(i) )/q(θ(i) ) ˜ IN := i=1 (biased). PN (i) (i) i=1 π(θ )/q(θ ) I One needs to know a good, heavy-tailed q. I Adaptive strategies for tuning q are possible [WKB+ 09]. R´ emi Bardenet Bayesian numerical methods First sampling methods 21 / 51 Importance Sampling 2/2 I If |h| ≤ M, then M2 Var(IˆN ) ≤ N I Z (π − q)2 dθ + 1 , q thus q has to be close to π and have heavier tails than π. Further: the q achieving the smallest variance is q=R |h|π . |h|πdθ I Better and achievable: For smaller variance, or if π is known up to a constant, use PN h(θ(i) )π(θ(i) )/q(θ(i) ) ˜ IN := i=1 (biased). PN (i) (i) i=1 π(θ )/q(θ ) I One needs to know a good, heavy-tailed q. I Adaptive strategies for tuning q are possible [WKB+ 09]. R´ emi Bardenet Bayesian numerical methods First sampling methods 21 / 51 Importance Sampling 2/2 I If |h| ≤ M, then M2 Var(IˆN ) ≤ N I Z (π − q)2 dθ + 1 , q thus q has to be close to π and have heavier tails than π. Further: the q achieving the smallest variance is q=R |h|π . |h|πdθ I Better and achievable: For smaller variance, or if π is known up to a constant, use PN h(θ(i) )π(θ(i) )/q(θ(i) ) ˜ IN := i=1 (biased). PN (i) (i) i=1 π(θ )/q(θ ) I One needs to know a good, heavy-tailed q. I Adaptive strategies for tuning q are possible [WKB+ 09]. R´ emi Bardenet Bayesian numerical methods First sampling methods 21 / 51 Going to high dimensions Benchmark Let π(θ) = N (0, Id ) and q(θ) = N (0, σId ). I Rejection sampling I I I needs σ ≥ 1. Fraction of accepted proposals goes as σ −d . Importance sampling I I √ yields infinite variance when σ ≤ 1/ 2, variance of the weights goes as d/2 σ2 . 2 2 − 1/σ R´ emi Bardenet Bayesian numerical methods First sampling methods 22 / 51 Going to high dimensions Benchmark Let π(θ) = N (0, Id ) and q(θ) = N (0, σId ). I Rejection sampling I I I needs σ ≥ 1. Fraction of accepted proposals goes as σ −d . Importance sampling I I √ yields infinite variance when σ ≤ 1/ 2, variance of the weights goes as d/2 σ2 . 2 2 − 1/σ R´ emi Bardenet Bayesian numerical methods First sampling methods 22 / 51 Going to high dimensions Benchmark Let π(θ) = N (0, Id ) and q(θ) = N (0, σId ). I Rejection sampling I I I needs σ ≥ 1. Fraction of accepted proposals goes as σ −d . Importance sampling I I √ yields infinite variance when σ ≤ 1/ 2, variance of the weights goes as d/2 σ2 . 2 2 − 1/σ R´ emi Bardenet Bayesian numerical methods First sampling methods 22 / 51 So far, so good ? We have seen that I Sums & integrals are ubiquitous in Statistics. I Numerical integration is limited to small dimensions. I “Basic” sampling is limited to full-information easy cases. I RS and IS, well-tuned, are efficient and versatile, Now what do we do when I it is not possible to sample from π directly, but only evaluate it pointwise, possibly up to a multiplicative constant: π(θ) = R p(x|θ)p(θ) . Θ p(x|θ)p(θ)dθ I We don’t know a good approximation q of π. I Θ is high-dimensional. Well, MCMC is bringing both answers and new problems ! R´ emi Bardenet Bayesian numerical methods First sampling methods 23 / 51 So far, so good ? We have seen that I Sums & integrals are ubiquitous in Statistics. I Numerical integration is limited to small dimensions. I “Basic” sampling is limited to full-information easy cases. I RS and IS, well-tuned, are efficient and versatile, Now what do we do when I it is not possible to sample from π directly, but only evaluate it pointwise, possibly up to a multiplicative constant: π(θ) = R p(x|θ)p(θ) . Θ p(x|θ)p(θ)dθ I We don’t know a good approximation q of π. I Θ is high-dimensional. Well, MCMC is bringing both answers and new problems ! R´ emi Bardenet Bayesian numerical methods First sampling methods 23 / 51 So far, so good ? We have seen that I Sums & integrals are ubiquitous in Statistics. I Numerical integration is limited to small dimensions. I “Basic” sampling is limited to full-information easy cases. I RS and IS, well-tuned, are efficient and versatile, Now what do we do when I it is not possible to sample from π directly, but only evaluate it pointwise, possibly up to a multiplicative constant: π(θ) = R p(x|θ)p(θ) . Θ p(x|θ)p(θ)dθ I We don’t know a good approximation q of π. I Θ is high-dimensional. Well, MCMC is bringing both answers and new problems ! R´ emi Bardenet Bayesian numerical methods First sampling methods 23 / 51 Summary 1 A very generic problem 2 First sampling methods 3 MCMC algorithms 4 A taste of a monster MCMC sampler for Auger R´ emi Bardenet Bayesian numerical methods MCMC algorithms 24 / 51 Metropolis’ algorithm 1/2 Goal is to explore Θ, spending more time in places where π is high. Figure: Taken from [Bis06] R´ emi Bardenet Bayesian numerical methods MCMC algorithms 25 / 51 Metropolis’ algorithm 2/2 I I Metropolis’ algorithm builds a random walk with symmetric steps q, that mimics independent draws from π. Symmetricity means q(x|y ) = q(y |x), e.g. q(y |x) = N (x, σ). MetropolisSampler π, q, T , θ(0) 1 S ← ∅. 2 for t ← 1 to T , 3 Sample θ∗ ∼ q(.|θ(t−1) ) and u ∼ U(0,1) . 4 Form the acceptance ratio ρ=1∧ π(θ∗ ) . π(θ(t−1) ) 5 if u < ρ, then θ(t) ← θ∗ else θ(t) ← θ(t−1) . 6 S ← S ∪ {θ(t) }. R´ emi Bardenet Bayesian numerical methods MCMC algorithms 26 / 51 The Metropolis-Hastings algorithm I When no symmetricity is assumed, change acceptance to π(y ) q(x|y ) . ρ(x, y ) = 1 ∧ q(y |x) π(x) I Remark the exploration/exploitation trade-off. MetropolisHastingsSampler π, q, T , θ(0) 1 S ← ∅. 2 for t ← 1 to T , 3 Sample θ∗ ∼ q(.|θ(t−1) ) and u ∼ U(0,1) . 4 Form the acceptance ratio ρ=1∧ π(θ∗ ) q(θ(t−1) |θ∗ ) . q(θ∗ |θ(t−1) ) π(θ(t−1) ) 5 if u < ρ, then θ(t) ← θ∗ else θ(t) ← θ(t−1) . 6 S ← S :: {θ(t) }. R´ emi Bardenet Bayesian numerical methods MCMC algorithms 27 / 51 All the maths in two slides 1/2 I A (stationary) Markov chain (θ(t) ) is defined through a kernel: P(θ(t+1) ∈ dy |θ(t) = x) = K (x, dy ). I If θ(0) = x, then P(θ (t) ∈ dy |θ (0) Z Z = x) = K (x, dx1 ) . . . K (dxt−1 , dy ) =: K t (x, dy ). R´ emi Bardenet Bayesian numerical methods MCMC algorithms 28 / 51 All the maths in two slides 2/2 I Under technical conditions, the MH kernel K satisfies kπ − K t (x, .)k → 0, ∀x. I This justifies a burn-in period after which samples are discarded. I Further results like a Law of Large Numbers guarantee that T 1 X h(θ(t) ) → T +1 Z h(θ)π(θ)d(θ) t=0 for h bounded. I Central Limit Theorem-type results also exist, see Section 6.7 of [RC04]. R´ emi Bardenet Bayesian numerical methods MCMC algorithms 29 / 51 Autocorrelation in equilibrium I For a scalar chain, define the integrated autocorrelation time X τint = 1 + 2 Corr (θ(t) , θ(t+k) ). k>0 I One can show that Var ! T 1 X τint h(θ(t) ) = Var h(θ(0) ) . T T t=1 Rule of thumb for the proposal Optimizing similar criteria leads to choosing q(.|x) = N (x, σ 2 ) s.t. I acceptance rate is ≈ 0.5 for d = 1, 2. I acceptance rate is ≈ 0.25 for d ≥ 3. R´ emi Bardenet Bayesian numerical methods MCMC algorithms 30 / 51 Autocorrelation in equilibrium I For a scalar chain, define the integrated autocorrelation time X τint = 1 + 2 Corr (θ(t) , θ(t+k) ). k>0 I One can show that Var ! T 1 X τint h(θ(t) ) = Var h(θ(0) ) . T T t=1 Rule of thumb for the proposal Optimizing similar criteria leads to choosing q(.|x) = N (x, σ 2 ) s.t. I acceptance rate is ≈ 0.5 for d = 1, 2. I acceptance rate is ≈ 0.25 for d ≥ 3. R´ emi Bardenet Bayesian numerical methods MCMC algorithms 30 / 51 Autocorrelation in equilibrium I For a scalar chain, define the integrated autocorrelation time X τint = 1 + 2 Corr (θ(t) , θ(t+k) ). k>0 I One can show that Var ! T 1 X τint h(θ(t) ) = Var h(θ(0) ) . T T t=1 Rule of thumb for the proposal Optimizing similar criteria leads to choosing q(.|x) = N (x, σ 2 ) s.t. I acceptance rate is ≈ 0.5 for d = 1, 2. I acceptance rate is ≈ 0.25 for d ≥ 3. R´ emi Bardenet Bayesian numerical methods MCMC algorithms 30 / 51 A practical wrap-up We have justified the acceptance ratio, the burn-in period, and the optimization of the proposal. MetropolisHastingsSampler π, q, T , θ(0) , Σ0 1 S ← ∅. 2 for t ← 1 to T , 3 Sample θ∗ ∼ N (.|θ(t−1) , σΣ0 ) and u ∼ U(0,1) . 4 Form the acceptance ratio ρ=1∧ π(θ∗ ) q(θ(t−1) |θ∗ ) . ∗ (t−1) q(θ |θ ) π(θ(t−1) ) 5 if u < ρ, then θ(t) ← θ∗ else θ(t) ← θ(t−1) . 6 if t ≤ Tb , 7 8 R´ emi Bardenet σ←σ+ 1 (acc. t 0.6 rate − 0.25/0.50). else if t > Tb , then S ← S ∪ θ(t) . Bayesian numerical methods MCMC algorithms 31 / 51 No maths on this slide Feel free to experiment with Laird Breyer’s applet on http://www.lbreyer.com/classic.html. R´ emi Bardenet Bayesian numerical methods MCMC algorithms 32 / 51 A case study from particle physics I Consider again our muonic signal reconstruction task, with p(data|θ) = N bins Y i=1 I n¯i (θ)ni −¯ni (θ) e . ni ! The model (the physics) suggest using specific independent priors for Aµ and tµ . 0.005 p 0.004 0.003 0.002 0.001 0.000 0 100 200 300 400 @ðPED 500 600 Figure: Prior on Aµ for a given zenith angle R´ emi Bardenet Bayesian numerical methods MCMC algorithms 33 / 51 Case study: simulation results log(A) chain 4 log(A) chain 4 marginal chain true value marginal chain true value 3.5 3.5 3 3 2.5 2.5 2 2 1.5 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 1.5 0 1000 2000 3000 4000 I One can often detect bad mixing by eye. I Acceptance for the left panel chain is only 3%. R´ emi Bardenet Bayesian numerical methods 5000 MCMC algorithms 6000 7000 8000 9000 10000 34 / 51 Case study: simulation results Marginal of logA 1500 Marginal of logA 2000 marginal histogram true value marginal mean marginal histogram true value marginal mean 1800 1600 1400 1000 1200 1000 800 500 600 400 200 0 1.5 2 2.5 3 3.5 4 0 1.5 2 2.5 I Marginals are simply component histograms ! I Even the marginals look ugly when mixing is bad. R´ emi Bardenet Bayesian numerical methods MCMC algorithms 3 3.5 4 34 / 51 Case study: simulation results Marginal of logA 1500 Marginal of logt 2500 marginal histogram true value marginal mean marginal histogram true value marginal mean 2000 1000 1500 1000 500 500 0 1.5 2 2.5 3 3.5 4 0 2.5 3 3.5 4 4.5 I From now on, we only consider the good mixing case. I Try to reproduce your marginals with different starting values! I Producing the chain was the hard part. Now everything is easy: estimation, credible intervals, ... R´ emi Bardenet Bayesian numerical methods MCMC algorithms 5 34 / 51 Summary 1 A very generic problem 2 First sampling methods 3 MCMC algorithms 4 A taste of a monster MCMC sampler for Auger A generative model for the tank signals [BK12] Reconstruction/Inference R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 35 / 51 Let’s go for Nµ > 1 I Recall Z ti pτ,td (t − tµ )dt. n¯i (Aµ , tµ ) = Aµ ti−1 50 45 40 35 ðPE 30 25 20 15 10 5 0 0 100 200 300 400 500 600 t @nsD I Now ni Poisson with mean n¯i (Aµ , tµ ) = Nµ X n¯i (Aµ j , tµ j ), j=1 R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 36 / 51 MCMC 101: Metropolis I Bayesian inference: obtain π(Aµ , tµ ) = p(Aµ , tµ |signal) ∝ p(signal|Aµ , tµ )p(Aµ , tµ ). I MCMC methods sample from π. R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 37 / 51 Metropolis (B) and adaptive Metropolis (B+G) algorithms ` ´ M π(x), X0 , Σ0 , T , 1 S ← ∅, Σ ← Σ0 2 for t ← 1 to T 3 ` ´ e ∼ N · |Xt−1 , Σ X 4 5 if e) π(X π(Xt−1 ) > U [0, 1] then Xt ← X 6 7 . proposal . accept else Xt ← Xt−1 8 S ← S ∪ {Xt } 9 . reject . update posterior sample 10 11 12 13 R´ emi Bardenet return S Bayesian numerical methods A taste of a monster MCMC sampler for Auger 38 / 51 Metropolis (B) and adaptive Metropolis (B+G) algorithms ` ´ AM π(x), X0 , Σ0 , T , µ0 , c 1 S ← ∅, 2 for t ← 1 to T 4 Σ ← cΣt−1 . scaled adaptive covariance ` ´ e X ∼ N · |Xt−1 , Σ . proposal 5 if 3 e) π(X π(Xt−1 ) > U [0, 1] then Xt ← X 6 7 Xt ← Xt−1 8 10 µt ← µt−1 + 11 Σt ← Σt−1 + 1` t . update posterior sample Xt − µt−1 1 `` t ´ Xt − µt−1 . update running mean and covariance ´` Xt − µt−1 ´| − Σt−1 ´ 1 (acc. rate − 0.25/0.50). c ← c + 0.6 t 12 R´ emi Bardenet . reject S ← S ∪ {Xt } 9 13 . accept else return S Bayesian numerical methods A taste of a monster MCMC sampler for Auger 38 / 51 Three problems, many “solutions” I Possibly high dimensions but also highly correlated model. I I I I I Use adaptive proposals. The number of muons Nµ is unknown. Use a nonparametric prior [Nea00] or use a Reversible Jump sampler [Gre95]. Likelihood P (n|Aµ , tµ ) is permutation invariant. I Indeed, if Nµ = 2, p(n| A1 , A2 , t1 , t2 ) = p(n| A2 , A1 , t2 , t1 ). I R´ emi Bardenet Marginals are useless, a problem known as label-switching [Ste00]. Bayesian numerical methods A taste of a monster MCMC sampler for Auger 39 / 51 Three problems, many “solutions” I Possibly high dimensions but also highly correlated model. I I I I I Use adaptive proposals. The number of muons Nµ is unknown. Use a nonparametric prior [Nea00] or use a Reversible Jump sampler [Gre95]. Likelihood P (n|Aµ , tµ ) is permutation invariant. I Indeed, if Nµ = 2, p(n| A1 , A2 , t1 , t2 ) = p(n| A2 , A1 , t2 , t1 ). I R´ emi Bardenet Marginals are useless, a problem known as label-switching [Ste00]. Bayesian numerical methods A taste of a monster MCMC sampler for Auger 39 / 51 Three problems, many “solutions” I Possibly high dimensions but also highly correlated model. I I I I I Use adaptive proposals. The number of muons Nµ is unknown. Use a nonparametric prior [Nea00] or use a Reversible Jump sampler [Gre95]. Likelihood P (n|Aµ , tµ ) is permutation invariant. I Indeed, if Nµ = 2, p(n| A1 , A2 , t1 , t2 ) = p(n| A2 , A1 , t2 , t1 ). I R´ emi Bardenet Marginals are useless, a problem known as label-switching [Ste00]. Bayesian numerical methods A taste of a monster MCMC sampler for Auger 39 / 51 Label-Switching I If the prior is also permutation invariant, then so is π(Aµ , tµ ). 300 250 200 150 100 50 0 0 R´ emi Bardenet 5000 10 000 15 000 Bayesian numerical methods 20 000 25 000 30 000 A taste of a monster MCMC sampler for Auger 40 / 51 Label-Switching I If the prior is also permutation invariant, then so is π(Aµ , tµ ). 50 45 40 35 ðPE 30 25 20 15 10 5 0 100 200 300 400 500 600 t @nsD R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 40 / 51 Label-Switching I If the prior is also permutation invariant, then so is π(Aµ , tµ ). 300 250 200 150 100 50 0 0 R´ emi Bardenet 5000 10 000 15 000 Bayesian numerical methods 20 000 25 000 30 000 A taste of a monster MCMC sampler for Auger 40 / 51 ` ´ AMOR π(x), X0 , T , µ0 , Σ0 , c 1 S ←∅ 2 for t ← 1 to T 5 Σ ← cΣt−1 . scaled adaptive covariance ` ´ e ∼ N · |Xt−1 , Σ X . proposal ` ´ e ∼ arg min L e P P X . pick an optimal permutation (µ ,Σ ) 6 e ←P eX e X 3 4 t−1 P∈P t−1 . permute ` ´ N PXt−1 |X , Σ P∈P if ` ´ > U [0, 1] then P π(Xt−1 ) P∈P N PX |Xt−1 , Σ e) π(X 7 P Xt ← X 8 9 . accept else Xt ← Xt−1 10 . reject 11 S ← S ∪ {Xt } 12 µt ← µt−1 + 13 Σt ← Σt−1 + 14 1 (acc. rate − 0.25/0.50). c ← c + 0.6 t 15 R´ emi Bardenet 1` t . update posterior sample Xt − µt−1 1 `` t ´ Xt − µt−1 . update running mean and covariance ´` Xt − µt−1 ´| − Σt−1 ´ return S Bayesian numerical methods A taste of a monster MCMC sampler for Auger 41 / 51 AMOR results on the previous example 50 45 40 35 ðPE 30 25 20 15 10 5 0 100 200 300 400 500 600 t @nsD R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 42 / 51 The integration methods ladder (adapted from Iain Murray’s) Growing item number means higher complexity and either higher efficiency or wider applicability. Check [RC04] when no further indication is given. 1 Quadrature, R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 43 / 51 The integration methods ladder (adapted from Iain Murray’s) Growing item number means higher complexity and either higher efficiency or wider applicability. Check [RC04] when no further indication is given. 1 Quadrature, 2 Quasi-MC, R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 43 / 51 The integration methods ladder (adapted from Iain Murray’s) Growing item number means higher complexity and either higher efficiency or wider applicability. Check [RC04] when no further indication is given. 1 Quadrature, 2 Quasi-MC, 3 Simple MC, Importance Sampling, R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 43 / 51 The integration methods ladder (adapted from Iain Murray’s) Growing item number means higher complexity and either higher efficiency or wider applicability. Check [RC04] when no further indication is given. 1 Quadrature, 2 Quasi-MC, 3 Simple MC, Importance Sampling, 4 MCMC (MH, Slice sampling, etc.), R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 43 / 51 The integration methods ladder (adapted from Iain Murray’s) Growing item number means higher complexity and either higher efficiency or wider applicability. Check [RC04] when no further indication is given. 1 Quadrature, 2 Quasi-MC, 3 Simple MC, Importance Sampling, 4 MCMC (MH, Slice sampling, etc.), 5 Adaptive MCMC [AFMP11], Hybrid MC [Nea10], Tempering methods [GRS96], SMC [DdFG01], particle MCMC [ADH10]. R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 43 / 51 The integration methods ladder (adapted from Iain Murray’s) Growing item number means higher complexity and either higher efficiency or wider applicability. Check [RC04] when no further indication is given. 1 Quadrature, 2 Quasi-MC, 3 Simple MC, Importance Sampling, 4 MCMC (MH, Slice sampling, etc.), 5 Adaptive MCMC [AFMP11], Hybrid MC [Nea10], Tempering methods [GRS96], SMC [DdFG01], particle MCMC [ADH10]. 6 ABC [FP12]. R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 43 / 51 Conclusions Take-home message I MC provides generic integration methods, I Potential applications in Physics are numerous: I I I I Producing a good mixing MCMC chain can be difficult Higher efficiency can result from: I I I in forward sampling (aka “simulation”), in Bayesian inference tasks. Learning dependencies. Exploiting existing/added structure of the problem. Broad range of methods allows to find the level of sophistication required by the difficulty of your problem. R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 44 / 51 Conclusions Take-home message I MC provides generic integration methods, I Potential applications in Physics are numerous: I I I I Producing a good mixing MCMC chain can be difficult Higher efficiency can result from: I I I in forward sampling (aka “simulation”), in Bayesian inference tasks. Learning dependencies. Exploiting existing/added structure of the problem. Broad range of methods allows to find the level of sophistication required by the difficulty of your problem. R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 44 / 51 “Introductory” references I Tutorials: Lots on videoconference.net. I I I. Murray’s video lessons: videolectures.net/mlss09uk murray mcmc/ A. Sokal’s tutorial, MCMC from a physicist’s point of view + applications to Statistical Physics [Sok96] I The Monte Carlo bible: C. Robert & G. Casella’s “Monte Carlo Methods” [RC04], and references within. I For more precise informations, please bug me. R´ emi Bardenet Bayesian numerical methods A taste of a monster MCMC sampler for Auger 45 / 51 Glossary In order of appearance: MAP MEP MC RS IS MCMC MH AMOR SMC ABC R´ emi Bardenet Maximum posterior (estimate) Mean posterior (estimate) Monte Carlo Rejection sampling Importance sampling Markov chain Monte Carlo Metropolis-Hastings (algorithm) Adaptive Metropolis with online relabeling Sequential Monte Carlo Approximate Bayesian computation Bayesian numerical methods A taste of a monster MCMC sampler for Auger 46 / 51 References I Christophe Andrieu, Arnaud Doucet, and Roman Holenstein, Particle Markov chain Monte Carlo methods, Journal of the Royal Statistical Society B (2010). Y. Atchad´e, G. Fort, E. Moulines, and P. Priouret, Bayesian time series models, ch. Adaptive Markov chain Monte Carlo : Theory and Methods, pp. 33–53, Cambridge Univ. Press, 2011. R. Bardenet, O. Capp´e, G. Fort, and B. K´egl, An adaptive Metropolis algorithm with online relabeling, International Conference on Artificial Intelligence and Statistics (AISTATS), (JMLR workshop and conference proceedings), vol. 22, April 2012, pp. 91–99. C. M. Bishop, Pattern recognition and machine learning, Springer, 2006. R´ emi Bardenet Bayesian numerical methods References 47 / 51 References II R. Bardenet and B. K´egl, An adaptive monte-carlo markov chain algorithm for inference from mixture signals, Proceedings of ACAT’11, Journal of Physics: Conference series, 2012. A. Doucet, N. de Freitas, and N. Gordon, Sequential Monte Carlo in practice, Springer, 2001. P. Fearnhead and D. Prangle, Constructing summary statistics for approx- imate bayesian computation; semi-automatic abc, Journal of the Royal Statistical Society B (2012). P. J. Green, Reversible jump Markov chain Monte Carlo computation and Bayesian model determination, Biometrika 82 (1995), no. 4, 711–732. W.R. Gilks, S. Richardson, and D. Spiegelhalter (eds.), Markov chain Monte Carlo in practice, Chapman & Hall, 1996. R´ emi Bardenet Bayesian numerical methods References 48 / 51 References III R. M. Neal, Markov chain sampling methods for Dirichlet process mixture models, Journal of Computational and Graphical Statistics 9 (2000), 249–265. , Handbook of Markov Chain Monte Carlo, ch. MCMC using Hamiltonian dynamics, Chapman & Hall / CRC Press, 2010. C. P. Robert and G. Casella, Monte Carlo statistical methods, Springer-Verlag, New York, 2004. A.D. Sokal, Monte Carlo methods in statistical mechanics: Foundations and new algorithms, 1996, Lecture notes at the Carg`ese summer school. M. Stephens, Dealing with label switching in mixture models, Journal of the Royal Statistical Society, Series B 62 (2000), 795–809. R´ emi Bardenet Bayesian numerical methods References 49 / 51 References IV D. Wraith, M. Kilbinger, K. Benabed, O. Capp´e, J.-F. Cardoso, G. Fort, S. Prunet, and C. P. Robert, Estimation of cosmological parameters using adaptive importance sampling, Phys. Rev. D 80 (2009), no. 2. R´ emi Bardenet Bayesian numerical methods References 50 / 51 Thanks for your attention R´ emi Bardenet Bayesian numerical methods References 51 / 51

© Copyright 2020