Bayesian parameter estimation in predictive engineering Damon - PowerPoint PPT Presentation

Bayesian parameter estimation in predictive engineering Damon McDougall Institute for Computational Engineering and Sciences, UT Austin 14th August 2014 1/27

Motivation Understand physical phenomena Observations of phenomena Mathematical model of phenomena (includes some parameters that characterise behaviour) Numerical model approximating mathematical model Find parameters in a situation of interest Use the parameters to do something cool 2/27

Understanding errors Reality errors Validation Mathematical model errors Verification Numerical model 3/27

Setup Model (usually a PDE): G ( u , θ ) where u is the initial condition and θ are model paramaters. u : perhaps an initial condition θ : perhaps some interesting model parameters (diffusion, convection speed, permeability field, material properties) Observations: i.i.d ∼ N (0 , σ 2 ) y j , k = u ( x j , t k ) + η j , k , η j , k η ∼ N (0 , σ 2 I ) y = G ( θ ) + η, ❀ Want: P ( θ | y ) ∝ P ( y | θ ) P ( θ ) Why? 4/27

Do we need Bayes’ theorem? Is Bayes’ theorem really necessary? We could minimise 2 σ 2 �G ( θ ) − y � 2 + 1 1 2 λ 2 � θ � 2 J ( θ ) = to get θ ∗ = argmin θ J ( θ ) 5/27

Do we need Bayes’ theorem? 6/27

Do we need Bayes’ theorem? 7/27

Do we need Bayes’ theorem? Bayesian methods involve estimating uncertainty (as well as mean). They’re equivalent. Deterministic optimisation: 1 1 2 σ 2 �G ( θ ) − y � 2 2 λ 2 � θ � 2 J ( θ ) = + � �� � � �� � misfit regularisation Bayesian framework: � � � � − 1 − 1 2 σ 2 �G ( θ ) − y � 2 2 λ 2 � θ � 2 exp( − J ( θ )) = exp exp � �� � � �� � likelihood prior = P ( y | θ ) P ( θ ) ∝ P ( θ | y ) 8/27

Method for solving Bayesian inverse problems • Kalman filtering/smoothing methods • Kalman filter (Kalman) • Ensemble Kalman filter (Evensen) • Variational methods • 3D VAR (Lorenc) • 4D VAR (Courtier, Talagrand, Lawless) • Particle methods • Particle filter (Doucet) • Sampling methods • Markov chain Monte Carlo (Metropolis, Hastings) This list is not exhaustive. The body of work is prodigious. 9/27

QUESO Nutshell: QUESO gives samples from P ( θ | y ) (called MCMC) • Library for Quantifying Uncertainty in Estimation, Simulation and Optimisation • Born in 2008 as part of PECOS PSAAP programme • Provides robust and scalable sampling algorithms for UQ in computational models • Open source • C++ • MPI for communication • Parallel chains, each chain can house several processes • Dependencies are MPI , Boost and GSL . Other optional features exist • https://github.com/libqueso/queso 10/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? 11/27

What does MCMC look like? � N E ( θ | y ) ≈ 1 k =1 θ k N 11/27

How to do MCMC? Sampling P ( θ | y ) i.i.d • Idea: Construct { θ k } ∞ k =1 cleverly such that { θ k } ∞ ∼ P ( θ | y ) k =1 1. Let θ j be the ‘current’ state in the sequence and construct a proposal , z 1 z = (1 − β 2 ) 2 θ j + βξ, some β ∈ (0 , 1) 2 σ 2 �G ( · ) − y � 2 1 2. Define Φ( · ) := 3. Compute α ( θ j , z ) = 1 ∧ exp(Φ( θ j ) − Φ( z )) 4. Let � θ with probability α ( θ j , z ) θ j +1 = with probability 1 − α ( θ j , z ) θ j • Take θ 1 to be a draw from P ( θ ) 12/27

How to do MCMC? Sampling P ( θ | y ) i.i.d • Idea: Construct { θ k } ∞ k =1 cleverly such that { θ k } ∞ ∼ P ( θ | y ) k =1 1. Let θ j be the ‘current’ state in the sequence and construct a proposal , z 1 z = (1 − β 2 ) 2 θ j + βξ, some β ∈ (0 , 1) 2 σ 2 �G ( · ) − y � 2 1 2. Define Φ( · ) := 3. Compute α ( θ j , z ) = 1 ∧ exp(Φ( θ j ) − Φ( z )) 4. Let � θ with probability α ( θ j , z ) θ j +1 = with probability 1 − α ( θ j , z ) θ j • Take θ 1 to be a draw from P ( θ ) 12/27

How to do MCMC? Sampling P ( θ | y ) i.i.d • Idea: Construct { θ k } ∞ k =1 cleverly such that { θ k } ∞ ∼ P ( θ | y ) k =1 1. Let θ j be the ‘current’ state in the sequence and construct a proposal , z 1 z = (1 − β 2 ) 2 θ j + βξ, some β ∈ (0 , 1) 2 σ 2 �G ( · ) − y � 2 1 2. Define Φ( · ) := 3. Compute α ( θ j , z ) = 1 ∧ exp(Φ( θ j ) − Φ( z )) 4. Let � θ with probability α ( θ j , z ) θ j +1 = with probability 1 − α ( θ j , z ) θ j • Take θ 1 to be a draw from P ( θ ) 12/27

How to do MCMC? Sampling P ( θ | y ) i.i.d • Idea: Construct { θ k } ∞ k =1 cleverly such that { θ k } ∞ ∼ P ( θ | y ) k =1 1. Let θ j be the ‘current’ state in the sequence and construct a proposal , z 1 z = (1 − β 2 ) 2 θ j + βξ, some β ∈ (0 , 1) 2 σ 2 �G ( · ) − y � 2 1 2. Define Φ( · ) := 3. Compute α ( θ j , z ) = 1 ∧ exp(Φ( θ j ) − Φ( z )) 4. Let � θ with probability α ( θ j , z ) θ j +1 = with probability 1 − α ( θ j , z ) θ j • Take θ 1 to be a draw from P ( θ ) 12/27

How to do MCMC? Sampling P ( θ | y ) i.i.d • Idea: Construct { θ k } ∞ k =1 cleverly such that { θ k } ∞ ∼ P ( θ | y ) k =1 1. Let θ j be the ‘current’ state in the sequence. Make a draw ξ ∼ P ( θ ) and construct a proposal , z 1 z = (1 − β 2 ) 2 θ j + βξ, some β ∈ (0 , 1) 2 σ 2 �G ( · ) − y � 2 1 2. Define Φ( · ) := 3. Compute α ( θ j , z ) = 1 ∧ exp(Φ( θ j ) − Φ( z )) 4. Let � θ with probability α ( θ j , z ) θ j +1 = with probability 1 − α ( θ j , z ) θ j • Take θ 1 to be a draw from P ( θ ) 12/27

How to do MCMC? Sampling P ( θ | y ) i.i.d • Idea: Construct { θ k } ∞ k =1 cleverly such that { θ k } ∞ ∼ P ( θ | y ) k =1 1. Let θ j be the ‘current’ state in the sequence. Make a draw ξ ∼ P ( θ ) and construct a proposal , z 1 z = (1 − β 2 ) 2 θ j + βξ, some β ∈ (0 , 1) 2 σ 2 �G ( · ) − y � 2 1 2. Define Φ( · ) := 3. Compute α ( θ j , z ) = 1 ∧ exp(Φ( θ j ) − Φ( z )) 4. Let � θ with probability α ( θ j , z ) θ j +1 = with probability 1 − α ( θ j , z ) θ j • Take θ 1 to be a draw from P ( θ ) 12/27

Why use QUESO? Other solutions are available, e.g. R, PyMC, emcee, MICA, Stan. QUESO solves the same problem, but: • Has been designed to be used with large forward problems • Has been used successfully with 5000+ cores • Leverages parallel MCMC algorithms • Supports for finite and infinite dimensional problems Statistical Application QUESO Forward Code 13/27

Why use QUESO? Chain 1 Chain 2 Chain 3 Node 1 Node 2 Node 3 Node 4 Node 5 Node 6 MCMC MCMC MCMC Node 1 Node 2 Node 3 Node 4 Node 5 Node 6 14/27

Example 1: convection-diffusion We are given a convection-diffusion model ( uc ) x − ( ν c x ) x = s , x ∈ [0 , 1] , c (0) = c (1) = 0 . Functions of x are: u , c and s . Constants are: ν (viscosity). The unkown is c , typically concentration. The underlying convection velocity is u . The forward problem: Given u and s , find c . 15/27

Example 1: convection-diffusion We are also given observations � ( uc ) x − ( ν c x ) x = s , x ∈ [0 , 1] , model c (0) = c (1) = 0 . � i.i.d ∼ N (0 , σ 2 ) , y j = c ( x j ) + η j , η j observations η ∼ N (0 , σ 2 I ) . ❀ y = G ( u ) + η, The observations are of c . We wish to learn about u . We will use Bayes’s theorem: P ( u | y ) ∝ P ( y | u ) P ( u ) True u = 1 − cos(2 π x ) True s = 2 π (1 − cos(2 π x )) cos(2 π x ) + 2 π sin 2 (2 π x ) + 4 π 2 ν sin(2 π x ) 16/27

Example 1: convection-diffusion How do we know we are solving the right PDE ( G ) to begin with? 10 0 10 − 2 10 − 4 10 − 6 L 2 error H 1 error 10 − 8 order 1 order 2 10 − 10 2 1 2 3 2 5 2 7 2 9 2 11 2 13 2 15 Number of grid points Note: Use the MASA [1] library to verify your forward problem. [1] Malaya et al., MASA: a library for verification using manufactured and analytical solutions, Engineering with Computers (2012) 17/27

Example 1: convection-diffusion Recap Bayes’s theorem, P ( u | y ) ∝ P ( y | u ) P ( u ) . Remember, we don’t know u but have observations and model: η ∼ N (0 , σ 2 I ) . y = G ( u ) + η, We also need a prior on u P ( u ) = N (0 , ( − ∆) − α ) . Aim is to get information from the posterior. 18/27

Example 1: convection-diffusion 1 . 2 1 . 0 0 . 8 0 . 6 0 . 4 0 . 2 � u k � L 2 � k � 1 i =0 u i � L 2 k 0 . 0 0 200 400 600 800 1000 Hundreds of iterations ( k ) 19/27

Example 1: convection-diffusion 0 . 20 � k 1 i =0 ( u i − � u i � k ) 2 � L 2 � k − 1 0 . 15 0 . 10 0 . 05 0 . 00 0 200 400 600 800 1000 Hundreds of iterations ( k ) 20/27