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 DS 2019, May 22

1 / 18

Outline

1

An SPDE from paleoclimate reconstruction Stochastic energy balance model State space model representation

2

Bayesian joint state-parameter estimation Sampling the posterior: Particle MCMC Ill-posedness: regularized posterior

3

Numerical results Parameter estimation State estimation

2 / 18

An SPDE from paleoclimate reconstruction

Paleoclimate: reconstruct past climate temperature from proxy data the temperature: a spatio-temporal process

◮ physically laws: energy balance → SPDEs ◮ discretized: a high-D process with spatial correlation

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 / 18

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

• gθ(u)

+W(t, x)

θ = (θk): unknown parameters

◮ prior: a range of physical values ◮ gθ(u) has a stable fixed point

W(t, x): Gaussian noise,

◮ white-in-time Matern-in-space

Data: sparse noisy observations

4 / 18

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

5 / 18

Joint parameter-state estimation SEBM: Un = g(θ, Un−1) + Wn Observation data: yn = HUn + Vn Goal: Given y1:N, we would like to jointly estimate (θ, U1:N) Gaussian prior for θ 12 spatial nodes, 100 time steps

• 1
• 0.5

0.5 1 1/

• 1
• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1 2/ 12 9 11 7 10 8 Solution at time step n =10 4 5 6 3 2 1 0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 20 40 60 80 100 Time step n 0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 un Trajectories of all 12 nodes

6 / 18

Bayesian 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 / 18

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 / 18

Ill-posed inverse problem For the Gaussian prior p(θ), unphysical samples of posterior: systems blowing up

9 / 18

Ill-posed inverse problem For the Gaussian prior p(θ), unphysical samples of posterior: systems blowing up Parameter estimation is ill-posed: Singular Fisher infomation matrix → large oscillation in sample θ from Gibbs θ| U1:N

2 3 4 5

Data size (log10N)

1 2 3

Std of errors (log10) Std of errors of MLE from noisy observations

1 4

10 / 18

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

= λ

• log p(θ) + 1

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

• λ = 1: Standard posterior N→∞

− − − − →∼ likelihood1 λ = N: regularized posterior pλ(θ, u1:N|y1:N) ∝ p(θ) [p(y1:N|u1:N)p(u1:N|θ)]1/N

1Bernstein-von Mises theorem 11 / 18

Parameter estimation

27 28 29 30 31 32 33 0.5 1 Posterior True value Posterior mean Prior

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

1

0.5 1 1.5

• 6.2
• 6
• 5.8
• 5.6
• 5.4
• 5.2
• 5
• 4.8
• 4.6

4

5 10

posterior close to prior; Errors in 100 simulations

θ0 θ1 θ4 Posterior mean

• 0.44 ± 0.58

0.09 ± 0.42 0.11 ± 0.20 MAP

• 0.32 ± 0.61

0.02 ± 0.42 0.03 ± 0.21

12 / 18

State estimation

Observed node: noise filtered Unobserved node: large spread

13 / 18

State estimation

Observed node: noise filtered Unobserved node: large spread

time = 20 0.95 1 1.05 State 0.1 0.2 0.3 0.4 0.5 0.6 Probability

Posterior True Mean

time = 60 0.95 1 1.05 State 0.1 0.2 0.3 0.4 0.5 0.6 time = 100 0.95 1 1.05 State 0.1 0.2 0.3 0.4 0.5 0.6 time = 20 0.9 1 1.1 State 0.1 0.2 0.3 0.4 0.5 Probability

Posterior True Mean

time = 60 0.9 1 1.1 State 0.1 0.2 0.3 0.4 0.5 time = 100 0.9 1 1.1 State 0.1 0.2 0.3 0.4 0.5

14 / 18

Observing 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 due to the need of regularization

15 / 18

Summary and discussion

Bayesian approach to jointly estimate parameter-state a stochastic energy balance model sparse and noisy data Ill-posed parameter estimation problem

(The parameters are correlated on a lowD manifold)

Introduced a regularized posterior: Enabling state estimation Large uncertainty in parameter estimation due to ill-posedness

16 / 18

Open questions

• 1. Re-parametrization/ nonparametric to avoid ill-posedness?
• 2. How many nodes need to be observed (for large mesh)?

(theory of determining modes)

17 / 18

Open questions

• 1. Re-parametrization/ nonparametric to avoid ill-posedness?
• 2. How many nodes need to be observed (for large mesh)?

(theory of determining modes)

Thank you!

18 / 18