Particle Markov Chain Monte Carlo Methods in Marine Biogeochemistry - - PowerPoint PPT Presentation

Particle Markov Chain Monte Carlo Methods in Marine Biogeochemistry

Lawrence Murray and Emlyn Jones CSIRO Mathematics, Informatics and Statistics

Lawrence Murray PMCMC in marine biogeochemistry : 2 of 41

Outline

• Case studies in marine biogeochemistry and generic challenges.
• Conventional state-space models
• Collapsed state-space models (to deal with the absence of a

closed-form transition density)

• The particle marginal Metropolis-Hastings (PMMH) sampler in

this context.

• Results and implications.

Acknowledgements: Eddy Campbell (CSIRO), John Parslow (CSIRO), Nugzar Margvelashvili (CSIRO), Noel Cressie (Ohio State University).

Lawrence Murray PMCMC in marine biogeochemistry : 3 of 41

Marine biogeochemical model

N

(DIN)

P

(Phytoplankton)

Z

(Zooplankton)

D

(Detritus)

Sea Surface Base of Mixed Layer Phytoplankton Growth Zooplankton Grazing Zooplankton Mortality Messy Feeding Remineralisation Mixing Sinking

Lawrence Murray PMCMC in marine biogeochemistry : 4 of 41

Issues, issues, issues...

• These models tend to have long memory and don’t mix well.
• Available observations may be sparse.
• The model may be strongly nonlinear, ...
• ...it might even be chaotic.
• The transition density is unlikely to have a closed form.

Lawrence Murray PMCMC in marine biogeochemistry : 5 of 41

Conventional state-space model

Y2 YT ... U1 U2 Y1 XT UT X2 X0 X1 ... ...

Lawrence Murray PMCMC in marine biogeochemistry : 6 of 41

Conventional state-space model

Sampling with sequential Monte Carlo (auxiliary particle filter)... For particle i at time t, extend xi

t−1 by drawing xi t ∼ q(xt) and

weight with:

wi

t = p(yt|xi t)p(xi t|xi t−1)

q(xi

t)

wi

t−1.

Lawrence Murray PMCMC in marine biogeochemistry : 7 of 41

Conventional state-space model

Sampling with sequential Monte Carlo (auxiliary particle filter)... For particle i at time t, extend xi

t−1 by drawing xi t ∼ q(xt) and

weight with:

wi

t = p(yt|xi t)p(xi t|xi t−1)

q(xi

t)

wi

t−1.

We don’t have that!

Lawrence Murray PMCMC in marine biogeochemistry : 8 of 41

Conventional state-space model

Y2 YT ... U1 U2 Y1 XT UT X2 X0 X1 ... ...

Lawrence Murray PMCMC in marine biogeochemistry : 9 of 41

Collapsed state-space model

YT ... U1 U2 UT Y1 X0 Y2 ...

Lawrence Murray PMCMC in marine biogeochemistry : 10 of 41

Collapsed state-space model

Sampling with sequential Monte Carlo (auxiliary particle filter)... For particle i at time t, extend ui

1:t−1 by drawing ui t ∼ q(ut) and

weight with:

wi

t = p(yt|ui 1:t, xi 0)p(ui t|ui 1:t−1, xi 0)

q(ui

t)

wi

t−1.

Lawrence Murray PMCMC in marine biogeochemistry : 11 of 41

Collapsed state-space model

Sampling with sequential Monte Carlo (auxiliary particle filter)... For particle i at time t, extend ui

1:t−1 by drawing ui t ∼ q(ut) and

weight with:

wi

t = p(yt|ui 1:t, xi 0)p(ui t)

q(ui

t)

wi

t−1.

We do have that!

Lawrence Murray PMCMC in marine biogeochemistry : 12 of 41

What should q(·) be?

• PF0: bootstrap
• PF1: as PF0 plus one-step single-pilot lookahead and resample.
• MUPF0: UKF run offline, use time marginals pN(ut|y1:t) as

proposals.

• MUPF1: as MUPF0 plus one-step single-pilot lookahead and

resample.

• CUPF0: UKF conditioned on each particle, so use pN(ut|u1:t−1, y1:t).
• CUPF1: as CUPF0 plus UKF lookahead and resample.

(UKF = Unscented Kalman Filter)

Lawrence Murray PMCMC in marine biogeochemistry : 13 of 41

Particle marginal Metropolis-Hastings

The target (posterior) density is π(u1:T, x0, θ|y1:T), factorised as either:

π1(u1:T, x0|θ, y1:T)π2(θ|y1:T)

• r

π1(u1:T|x0, θ, y1:T)π2(x0, θ|y1:T).

In either case π1 is targeted with an auxiliary particle filter, π2 with Metropolis-Hastings [Andrieu et al., 2010, Pitt et al., 2011]. The second factorisation should be preferred when prior information

• ver x0 is scarce, so that importance sampling of it will be degenerate.

Lawrence Murray PMCMC in marine biogeochemistry : 14 of 41

PZ model

Simple phytoplankton-zooplankton model [Jones et al., 2010]. Lotka-Volterra with stochastic growth rate and quadratic mortality term.

dP dt = αtP − cPZ dZ dt = ecPZ − mlZ − mqZ2 αt ∼ N(µ, σ). P is observed with noise. Simulated data used here.

Lawrence Murray PMCMC in marine biogeochemistry : 15 of 41

PZ model: state estimate

1 2 3 4 5 6 7 20 40 60 80 100 P t Prior Posterior Observed Truth 0.5 1 1.5 2 2.5 20 40 60 80 100 Z

• 2
• 1

1 2 20 40 60 80 100 α t

Lawrence Murray PMCMC in marine biogeochemistry : 16 of 41

NPZD model

The model features nine noise terms, ξi for i = 1, . . . , 9, each coupled to a univariate autoregressive process Bi.

Bi(t + ∆t) = Bi(t) · (1 − ∆t/τP) + (µi + PDF · σiξi) · ∆t/τP,

where ∆t is a discrete time step (one day), µi a parameter to be estimated, PDF a common diversity factor to be estimated, σi a prescribed scaling factor, and τP a common characteristic time scale, also prescribed.

Lawrence Murray PMCMC in marine biogeochemistry : 17 of 41

NPZD model

A multiplicative temperature correction Tc is applied to all rate processes, for which a Q10 formulation for dependence on temperature, T , is used:

Tc = Q(T−Tref)/10

10

,

where Tref is a reference temperature, and Q10 a prescribed constant.

Lawrence Murray PMCMC in marine biogeochemistry : 18 of 41

NPZD model

The zooplankton grazing rate (gr) is dependent on the phytoplankton concentration (zooplankton functional response):

gr = Tc · IZ · Aυ (1 + Aυ) ,

(1) where υ is a given power.

Lawrence Murray PMCMC in marine biogeochemistry : 19 of 41

NPZD model

The relative availability of phytoplankton, A, is

A = ClZ · P IZ ; IZ is the maximum zooplankton ingestion rate (mg P per mg Z

per day); ClZ is the maximum clearance rate (volume in m3 swept clear per mg Z per day). For υ = 1, this takes the form of a Type-2 functional response (standard rectangular hyperbola) [Holling, 1966], and for υ > 1 a Type-3 sigmoid functional response.

Lawrence Murray PMCMC in marine biogeochemistry : 20 of 41

NPZD model

A quadratic formulation for zooplankton mortality is adopted after Steele [1976] and Steele and Henderson [1992]:

m = Tc · mQ · Z,

where the quadratic mortality rate mQ has units of d−1(mgZm−3)−1.

Lawrence Murray PMCMC in marine biogeochemistry : 21 of 41

NPZD model

The detrital remineralization rate is dependent only on temperature:

r = Tc · rD,

where rD prescribes the remineralisation rate at a reference temperature.

Lawrence Murray PMCMC in marine biogeochemistry : 22 of 41

NPZD model

The phytoplankton specific growth rate, g, depends on temperature,

T , available light or irradiance, E, and dissolved inorganic nutrient, N: g = Tc · gmax · hE · hN/(hE + hN).

Lawrence Murray PMCMC in marine biogeochemistry : 23 of 41

NPZD model

The light-limitation factor is given by

hE = 1 − exp(−α · λmax · E/gmax),

where α is the initial slope of the photosynthesis versus irradiance curve (mg C mg Chla−1 mol photon−1 m2), and λmax is the maximum

Chla : C ratio (mg Chla mg C−1).

Lawrence Murray PMCMC in marine biogeochemistry : 24 of 41

NPZD model

E is the mean photosynthetic available radiation (PAR) in the mixed

layer and is given by

E = E0.(1 − exp(−Kz))/Kz,

where E0 is the mean daily photosynthetically available radiation (PAR) just below the air-sea interface and Kz is given by

Kz = (KW + aCh · Chla) · MLD.

(2)

KW is attenuation due to the seawater and aCh.

Lawrence Murray PMCMC in marine biogeochemistry : 25 of 41

NPZD model

The nutrient-limitation factor is given by

hN = N (gmax · Tc/aN) + N ,

where aN is the maximum specific affinity for nitrogen uptake (d−1 mg N −1 m3).

Lawrence Murray PMCMC in marine biogeochemistry : 26 of 41

NPZD model

The phytoplankton N : C ratio, χ, predicted by the model is given by

χ = χmin · hE + χmax · hN hE + hN ,

where χmin and χmax are the prescribed minimum and maximum

N : C ratios (mg N mg C−1).

Lawrence Murray PMCMC in marine biogeochemistry : 27 of 41

NPZD model

The equations governing interactions between the remaining state variables {P, Z, D, N} are:

dP dt = g · P − gr · Z + κ MLD · (BCP − P) dZ dt = EZ · gr · Z − m · Z dD dt = (1 − EZ) · fD · gr · Z + m · Z − r · D − SD · D MLD + κ MLD · (BCD − D) dN dt = −g · P + (1 − EZ) · (1 − fD) · gr · Z + r · D + κ MLD · (BCN − N).

Lawrence Murray PMCMC in marine biogeochemistry : 28 of 41

State estimation (observed)

1971 1972 1973 1974 1975 50 100 150 200 250 300 350

DIN (µ g N l−1)

1971 1972 1973 1974 1975 −8 −6 −4 −2 2 4

ln(Chla (µ g Chla l−1))

Prior−95% Posterior−95% Observations Prior−median Posterior−median

Lawrence Murray PMCMC in marine biogeochemistry : 29 of 41

State estimation (unobserved)

1971 1972 1973 1974 1975 −4 −2 2 4 6

ln(P(µ g N l−1)

1971 1972 1973 1974 1975 −2 −1 1 2 3

ln(Z(µ g N l−1)

Prior−95% Posterior−95% Prior−median Posterior−median 1971 1972 1973 1974 1975 −6 −4 −2 2 4

ln(D(µ g N l−1)

Lawrence Murray PMCMC in marine biogeochemistry : 30 of 41

Parameter estimation

−4 −3.5 −3 1 2 ln(KW) −4.5 −4 −3.5 0.5 1 ln(aCh) 2 4 6 8 0.2 0.4 sD −0.8 −0.6 −0.4 2 4 ln(fD) −2.5 −2 −1.5 −1 −0.5 0.5 1 ln(PDF) −2.5 −2 −1.5 −1 −0.5 0.5 1 ln(ZDF) −1 1 2 0.2 0.4 0.6 0.8 ln(µg

max)

−4.5 −4 −3.5 −3 −2.5 0.5 1 ln(µλ

max)

−2 −1.5 −1 −0.5 0.5 1 ln(µR

N

) −4 −2 0.2 0.4 ln(µa

N

) 1 2 3 0.2 0.4 0.6 ln(µI

Z

) −4 −2 2 0.2 0.4 0.6 ln(µCl

Z

) −1.5 −1 −0.5 0.5 1 1.5 ln(µE

Z

) −3 −2 −1 0.2 0.4 0.6 0.8 ln(µr

D

) −6 −4 −2 0.2 0.4 ln(µm

Q

) Posterior Prior

Lawrence Murray PMCMC in marine biogeochemistry : 31 of 41

State forecast

07/74 02/75 09/75 50 100 150 200 250

DIN(µ g N l

−1)

07/74 02/75 09/75 −4 −2 2 4 6

ln(P(µ g N l

−1))

07/74 02/75 09/75 −2 −1 1 2 3 4

ln(Z(µ g N l

−1))

07/74 02/75 09/75 −8 −6 −4 −2 2 4 6

ln(D(µ g N l

−1))

07/74 02/75 09/75 −8 −6 −4 −2 2 4

ln(Chla(µ g N l

−1))

Prior−95% Posterior−95% Forecast−95% Observations Prior−median Posterior−median Forecast−median

Lawrence Murray PMCMC in marine biogeochemistry : 32 of 41

PZ model: convergence

1 1.002 1.004 1.006 1.008 1.01 1.012 5000 10000 15000 20000 Rp Step PF0 MUPF0 CUPF0 PF1 MUPF1 CUPF1 1 1.002 1.004 1.006 1.008 1.01 1.012 5000 10000 15000 20000 Step PF0 MUPF0 CUPF0 PF1 MUPF1 CUPF1

Convergence rates of Markov chains for the PZ case study, (left) particle-matched, and (right) compute-matched. Each line shows the evolution of the ˆ

Rp statistic of Brooks and Gelman

[1998] for a particular method.

Lawrence Murray PMCMC in marine biogeochemistry : 33 of 41

NPZD model: convergence

1 1.2 1.4 1.6 1.8 2 2.2 2.4 2000 4000 6000 8000 10000 Step PF0 MUPF0 CUPF0 PF1 MUPF1 CUPF1 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2000 4000 6000 8000 10000 Step PF0 MUPF0 CUPF0 PF1 MUPF1 CUPF1

Convergence rates of Markov chains for the NPZD case study, (left) particle-matched, and (right) compute-matched. Each line shows the evolution of the ˆ

Rp statistic of Brooks and

Gelman [1998] for a particular method.

Lawrence Murray PMCMC in marine biogeochemistry : 34 of 41

PZ model: acceptance rates

σ PF0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 MUPF0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 CUPF0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 σ µ PF1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 µ MUPF1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 µ CUPF1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

Lawrence Murray PMCMC in marine biogeochemistry : 35 of 41

PZ model: acceptance rates, compute-matched

σ PF0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 MUPF0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 CUPF0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 σ µ PF1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 µ MUPF1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 µ CUPF1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

Lawrence Murray PMCMC in marine biogeochemistry : 36 of 41

NPZD model: acceptance rates

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Cumulative density Conditional acceptance rate (CAR) PF0 MUPF0 CUPF0 PF1 MUPF1 CUPF1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 Cumulative density Conditional acceptance rate (CAR) PF0 MUPF0 CUPF0 PF1 MUPF1 CUPF1

Empirical cdfs of the acceptance rates for each method on the NPZD case study, (left) particle-matched, and (right) compute-matched.

Lawrence Murray PMCMC in marine biogeochemistry : 37 of 41

NPZD model: parameters vs acceptance rates

KW 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 aCh 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 SD 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 fD 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 PDF 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 ZDF 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µgmax 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µλmax 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µRN 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µaN 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µIZ 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µClZ 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µEZ 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µrD 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 µmQ 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

Lawrence Murray PMCMC in marine biogeochemistry : 38 of 41

Lessons for PMMH

• Initialisation is extremely important, sampler will fail if initialised

in a region where the acceptance rate is low.

• A particularly informative prior distribution may bias the posterior

into regions where the acceptance rate is low.

• Is the rule-of-thumb 23% acceptance rate appropriate in this

context?

Lawrence Murray PMCMC in marine biogeochemistry : 39 of 41

Summary

• Environmental models such as these in marine biogeochemistry

pose significant challenges to Monte Carlo methodology.

• One particular challenge is the absence of a closed-form transition
• density. The collapsed state-space model facilitates one means
• f improving sampler performance without needing this.
• The particle marginal Metropolis-Hastings sampler can have

significant mixing and acceptance rate problems, initialisation and careful consideration of the prior distribution are important.

Lawrence Murray PMCMC in marine biogeochemistry : 40 of 41

References

• C. Andrieu, A. Doucet, and R. Holenstein. Particle Markov chain Monte Carlo methods. Journal
• f the Royal Statistical Society Series B, 72:269–302, 2010.
• S. P

. Brooks and A. Gelman. General methods for monitoring convergence of iterative

• simulations. Journal of Computational and Graphical Statistics, 7:434–455, 1998.
• C. S. Holling. The functional response of predators to prey density and its role in mimicry and

population regulation. Memoirs of the Entomology Society of Canada, 45:60, 1966.

• E. Jones, J. Parslow, and L. M. Murray. A Bayesian approach to state and parameter estimation

in a phytoplankton-zooplankton model. Australian Meteorological and Oceanographic Journal, 59(SP):7–16, 2010.

• M. K. Pitt, R. S. Silva, P

. Giordani, and R. Kohn. Auxiliary particle filtering within adaptive Metropolis-Hastings sampling. 2011. URL http://arxiv.org/abs/1006.1914.

• J. Steele. Role of predation in ecosystem models. Marine Biology, 35(1):9–11, 1976. ISSN

0025-3162.

• J. Steele and E. Henderson. The role of predation in plankton models. Journal of Plankton

Research, 14(1):157–172, 1992.

CSIRO Mathematics, Informatics and Statistics Lawrence Murray Phone: +61 8 9333 6480 Email: lawrence.murray@csiro.au Web: www.cmis.csiro.au Contact Us Phone: 1300 363 400 or +61 3 9545 2176 Email: enquiries@csiro.au Web: www.csiro.au