Joint stateparameter estimation for nonlinear stochastic energy - - PowerPoint PPT Presentation

Joint state–parameter estimation for nonlinear stochastic energy balance models

Fei Lu 1 Nils Weitzel2 Adam Monahan3

1Department of Mathematics, Johns Hopkins 2Meteorological Institute, University of Bonn, Germany 3School of Earth and Ocean Science, University of Victoria, Canada

SIAM Minisymposium on Data Assimilation: Theory and Practice JMM 2019, January 17

1 / 20

Outline

1

Motivation Stochastic energy balance model State space model representation

2

Bayesian inference Particle MCMC Regularized posterior

3

Numerical study 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,

• cean 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)

∂tu = QT

• transport

+ QSW

• absorbed

shortwave

+ QSH

• sensible

heat

+ QLH

• latent

heat

+ QLW

• longwave

surf.→atmos.

− QLPW

longwave into space

= ∇ · (ν∇u) + θ0 + θ1u + θ4u4

• +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: yti =

• Ai

u(ti, x) dx + Vi, {Ai} are regions/locations of observations Gaussian noise {Vi}, iid, variance known Linear operator in state u

5 / 20

State space model formulation

SEBM: ∂tu = ∇ · (ν∇u) +

• k=0,1,4

θkuk + W(t, x) Observation data: yti = H(u(ti, x)) + Vi Discretization (simplification): finite elements in space semi-backward Euler in time State space model SEBM: Un = g(θ, Un−1) + Wn Observation data: Yn = HUn + Vn Goal: Given y1:N, we would like to jointly estimate (θ, U1:N) Bayesian approach to quantify uncertainty

6 / 20

Joint state-parameter estimation

Bayesian approach: p(θ, u1:N|y1:N) ∝ p(θ)p(u1:N|θ)p(y1:N|u1:N) Posterior: quantifies the uncertainties Approximate the posterior by sampling high dimensional (> 103 ), non-Gaussian, mixed types of variables θ, u1:N Gibbs Monte Carlo: U1:N|θ and θ|U iteration

◮ U1: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 θ| U1:N

10 / 20

Regularized posterior

Recall the regularization in variational approach Variational: ( θ, u1:N) = arg min

(θ,u1:n)

Cλ,y1:N(θ, u1:N) Bayesian : pλ(θ, u1:N|y1:N) ∝ p(θ)λp(y1:N|u1:N)p(u1:N|θ)

Cλ,y1:N(θ, u1:N) = λ log p(θ)

• regularization

+ log[p(y1:N|u1:N)p(u1:N|θ)]

• likelihood

λ = 1: Standard posterior N→∞ − − − − →∼ likelihood λ = N: regularized posterior

11 / 20

Numerical tests

Physical parameter set up: Gaussian prior

θ0 θ1 θ4 mean 30.11

• 24.08
• 5.40

std 0.82 0.46 0.20

temperature near an equilibrium point (normalized, ≈ 1)

20 40 60 80 100 120 0.96 0.98 1 1.02 1.04 1.06 All nodes 2 4 6 8 10 12 0.96 0.98 1 1.02 1.04 1.06 Mean and std of Utrue Mean std

Dimension of the states: 1200 12 spatial nodes 100 time steps

• bserve 6 nodes each time;

Randomly generate a true value for parameter from prior

12 / 20

Diagnosis of Markov Chain

Chain length: 1000 (with 30% burnin)

50 100 150 200 Time Lag

• 0.4
• 0.2

0.2 0.4 0.6 0.8 1 Correlation of chain: states 50 100 150 200 Time Lag

• 0.4
• 0.2

0.2 0.4 0.6 0.8 1 Correlation of chain: parameters

Correlation length: 10-30 steps

20 40 60 80 100 120 Time 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Update rates Update (Acceptance) rate of the MCMC

State update rate: > 0.55

13 / 20

Parameter estimation

24 25 26 27 28 29 30 31 32 33 34 0.5 1

Marginal posterior posterior true value posterior mean prior mean/std

• 25.5
• 25
• 24.5
• 24
• 23.5
• 23
• 22.5
• 22

0.5 1 1.5

Marginal posterior

1

• 5.8
• 5.6
• 5.4
• 5.2
• 5
• 4.8
• 4.6
• 4.4
• 4.2
• 4

2 4 6

Marginal posterior

4

Marginals of the posterior

θ0, θ1 OK, +bias θ4 posterior close to prior Errors in 100 simulations

θ0 θ1 θ4 Mean

• 0.74

0.11 0.22 Std 0.73 0.46 0.20

◮ -bias θ0, +bias in θ4 14 / 20

State estimation

Ensemble of sample trajectories:

Observed node: filtered out noise Unobserved node: large spread, mean close to truth

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