Joint state–parameter estimation for nonlinear stochastic energy balance models Fei Lu 1 Nils Weitzel 2 Adam Monahan 3 1 Department of Mathematics, Johns Hopkins 2 Meteorological Institute, University of Bonn, Germany 3 School of Earth and Ocean Science, University of Victoria, Canada SIAM Minisymposium on Data Assimilation: Theory and Practice JMM 2019, January 17 1 / 20
Outline Motivation 1 Stochastic energy balance model State space model representation Bayesian inference 2 Particle MCMC Regularized posterior Numerical study 3 Diagnosis of Markov Chain Parameter estimation State estimation 2 / 20
Motivation Paleoclimate : reconstruct past climate temperature from proxy data Spatio-temporal evolution ◮ spatial correlations ◮ physically laws: energy balance → SPDEs Sparse and noisy data ◮ Proxy data: historical data, tree rings, ice cores, fossil pollen, ocean sediments, coral etc. Plan: inference of SPDEs from sparse noisy data joint state-parameter estimation 3 / 20
The SPDEs: stochastic Energy Balance Models Idealized atmospheric energy balance ( Fanning&Weaver1996 ) ∂ t u = Q T + Q SW + Q SH + Q LH + Q LW − Q LPW ���� ���� ���� ���� ���� � �� � transport absorbed sensible latent longwave longwave shortwave heat heat surf. → atmos. into space ∇ · ( ν ∇ u ) + θ 0 + θ 1 u + θ 4 u 4 = � + W ( t , x ) � �� u ( t , x ) normalized temperature ( ≈ 1) θ = ( θ k ) : unknown parameters: ◮ prior: a range of physical values W ( t , x ) : Gaussian noise, ◮ white-in-time Matern-in-space 4 / 20
Data: observation model Observation at sparse locations/regions: � y t i = u ( t i , x ) dx + V i , A i { A i } are regions/locations of observations Gaussian noise { V i } , iid, variance known Linear operator in state u 5 / 20
State space model formulation � θ k u k + W ( t , x ) SEBM: ∂ t u = ∇ · ( ν ∇ u ) + k = 0 , 1 , 4 Observation data: y t i = H ( u ( t i , x )) + V i Discretization (simplification): finite elements in space semi-backward Euler in time State space model SEBM: U n = g ( θ, U n − 1 ) + W n Observation data: Y n = HU n + V n Goal: Given y 1 : N , we would like to jointly estimate ( θ, U 1 : N ) Bayesian approach to quantify uncertainty 6 / 20
Joint state-parameter estimation Bayesian approach: p ( θ, u 1 : N | y 1 : N ) ∝ p ( θ ) p ( u 1 : N | θ ) p ( y 1 : N | u 1 : N ) Posterior: quantifies the uncertainties Approximate the posterior by sampling high dimensional ( > 10 3 ), non-Gaussian, mixed types of variables θ, u 1 : N Gibbs Monte Carlo: U 1 : N | θ and θ | U iteration ◮ U 1 : N | θ needs highD proposal density → Sequential MC ◮ combine SMC with Gibbs (MCMC) → Particle MCMC methods based on conditional SMC 7 / 20
Sampling: particle MCMC Particle MCMC (Andrieu&Doucet&Holenstein10) Combines Sequential MC with MCMC: ◮ SMC: seq. importance sampling → highD proposal density ◮ conditional SMC: keep a reference trajectory in SMC ◮ MCMC transition by conditional SMC → target distr invariant even w/ a few particles Particle Gibbs with Ancestor Sampling (Lindsten&Jordan&Schon14) ◮ Update the ancestor of the reference trajectory ◮ Improving mixing of the chain 8 / 20
However, standard Bayesian approach does not work: for a Gaussian prior p ( θ ) , unphysical samples of posterior: systems blowing up 9 / 20
However, standard Bayesian approach does not work: for a Gaussian prior p ( θ ) , unphysical samples of posterior: systems blowing up Parameter estimation is ill-posed: Singular Fisher infomation matrix for full perfect observation → large oscillation in sample θ from Gibbs θ | � U 1 : N 10 / 20
Regularized posterior Recall the regularization in variational approach ( � θ, � Variational: u 1 : N ) = arg min C λ, y 1 : N ( θ, u 1 : N ) ( θ, u 1 : n ) p λ ( θ, u 1 : N | y 1 : N ) ∝ p ( θ ) λ p ( y 1 : N | u 1 : N ) p ( u 1 : N | θ ) Bayesian : C λ, y 1 : N ( θ, u 1 : N ) = λ log p ( θ ) + log[ p ( y 1 : N | u 1 : N ) p ( u 1 : N | θ )] � �� � � �� � regularization likelihood λ = 1: Standard posterior N →∞ − − − − →∼ likelihood λ = N : regularized posterior 11 / 20
Numerical tests Physical parameter set up: θ 0 θ 1 θ 4 Gaussian prior mean 30.11 -24.08 -5.40 std 0.82 0.46 0.20 temperature near an equilibrium point (normalized, ≈ 1) All nodes 1.06 1.04 Dimension of the states: 1200 1.02 12 spatial nodes 1 0.98 100 time steps 0.96 0 20 40 60 80 100 120 observe 6 nodes each time; Mean and std of Utrue 1.06 Mean std Randomly generate a true value 1.04 1.02 for parameter from prior 1 0.98 0.96 0 2 4 6 8 10 12 12 / 20
Diagnosis of Markov Chain Chain length: 1000 (with 30% burnin) Correlation of chain: states Correlation of chain: parameters 1 1 Update (Acceptance) rate of the MCMC 1 0.8 0.8 0.95 0.6 0.6 0.9 0.85 0.4 0.4 Update rates 0.8 0.75 0.2 0.2 0.7 0 0 0.65 0.6 -0.2 -0.2 0.55 0 20 40 60 80 100 120 Time -0.4 -0.4 0 50 100 150 200 0 50 100 150 200 Time Lag Time Lag State update rate: > 0.55 Correlation length: 10-30 steps 13 / 20
Parameter estimation Marginal posterior 0 1 posterior true value 0.5 posterior mean prior mean/std θ 0 , θ 1 OK, +bias θ 4 0 24 25 26 27 28 29 30 31 32 33 34 Marginal posterior 1 posterior close to prior 1.5 1 0.5 Errors in 100 simulations 0 θ 0 θ 1 θ 4 -25.5 -25 -24.5 -24 -23.5 -23 -22.5 -22 Mean -0.74 0.11 0.22 Marginal posterior 4 6 Std 0.73 0.46 0.20 4 ◮ -bias θ 0 , +bias in θ 4 2 0 -5.8 -5.6 -5.4 -5.2 -5 -4.8 -4.6 -4.4 -4.2 -4 Marginals of the posterior 14 / 20
State estimation Ensemble of sample trajectories: Unobserved node: Observed node: large spread, mean close to truth filtered out noise 15 / 20
Observe more or less nodes When more modes are observed: State estimation gets more accurate Parameter estimation does not improve much: the posterior keeps close to prior. 16 / 20
Summary and discussion Bayesian approach to jointly estimate parameter-state a stochastic energy balance model sparse and noisy data Estimate both parameters and states ◮ regularized posterior due to singular Fisher matrix ◮ Gibbs sampling via PGAS 17 / 20
Summary and discussion Bayesian approach to jointly estimate parameter-state a stochastic energy balance model sparse and noisy data Estimate both parameters and states ◮ regularized posterior due to singular Fisher matrix ◮ Gibbs sampling via PGAS Results: State estimation: ◮ filtered noise on observed nodes; ◮ large uncertainty in unobserved modes Parameter estimation: ◮ slightly biased estimators ◮ posterior close to prior 18 / 20
Open quetions 1. Re-parametrization to avoid singular Fish matrix? 2. How many nodes need to be observed (for large mesh)? (theory of determining modes) 19 / 20
Open quetions 1. Re-parametrization to avoid singular Fish matrix? 2. How many nodes need to be observed (for large mesh)? (theory of determining modes) Thank you! 20 / 20
Recommend
More recommend
