Particle Markov Chain Monte Carlo Methods in Marine Biogeochemistry Lawrence Murray and Emlyn Jones CSIRO Mathematics, Informatics and Statistics
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 : 2 of 41
Marine biogeochemical model Sea Surface P Z Zooplankton Grazing (Phytoplankton) (Zooplankton) Phytoplankton Zooplankton Mortality Growth Messy Feeding N D Remineralisation (DIN) (Detritus) Base of Mixed Layer Mixing Sinking Lawrence Murray PMCMC in marine biogeochemistry : 3 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 : 4 of 41
Conventional state-space model ... U 1 U 2 U T ... X 0 X 1 X 2 X T ... Y 1 Y 2 Y T Lawrence Murray PMCMC in marine biogeochemistry : 5 of 41
Conventional state-space model Sampling with sequential Monte Carlo (auxiliary particle filter)... For particle i at time t , extend x i t − 1 by drawing x i t ∼ q ( x t ) and weight with: t = p ( y t | x i t ) p ( x i t | x i t − 1 ) w i w i t − 1 . q ( x i t ) 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 x i t − 1 by drawing x i t ∼ q ( x t ) and weight with: t = p ( y t | x i t ) p ( x i t | x i t − 1 ) w i w i t − 1 . q ( x i t ) We don’t have that! Lawrence Murray PMCMC in marine biogeochemistry : 7 of 41
Conventional state-space model ... U 1 U 2 U T ... X 0 X 1 X 2 X T ... Y 1 Y 2 Y T Lawrence Murray PMCMC in marine biogeochemistry : 8 of 41
Collapsed state-space model ... U 1 U 2 U T X 0 ... Y 1 Y 2 Y T Lawrence Murray PMCMC in marine biogeochemistry : 9 of 41
Collapsed state-space model Sampling with sequential Monte Carlo (auxiliary particle filter)... For particle i at time t , extend u i 1: t − 1 by drawing u i t ∼ q ( u t ) and weight with: t = p ( y t | u i 1: t , x i 0 ) p ( u i t | u i 1: t − 1 , x i 0 ) w i w i t − 1 . q ( u i t ) 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 u i 1: t − 1 by drawing u i t ∼ q ( u t ) and weight with: t = p ( y t | u i 1: t , x i 0 ) p ( u i t ) w i w i t − 1 . q ( u i t ) We do have that! Lawrence Murray PMCMC in marine biogeochemistry : 11 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 p N ( u t | y 1: t ) as proposals. • MUPF1: as MUPF0 plus one-step single-pilot lookahead and resample. • CUPF0: UKF conditioned on each particle, so use p N ( u t | u 1: t − 1 , y 1: t ) . • CUPF1: as CUPF0 plus UKF lookahead and resample. (UKF = Unscented Kalman Filter) Lawrence Murray PMCMC in marine biogeochemistry : 12 of 41
Particle marginal Metropolis-Hastings The target (posterior) density is π ( u 1: T , x 0 , θ | y 1: T ) , factorised as either: π 1 ( u 1: T , x 0 | θ , y 1: T ) π 2 ( θ | y 1: T ) or π 1 ( u 1: T | x 0 , θ , y 1: T ) π 2 ( x 0 , θ | y 1: 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 over x 0 is scarce, so that importance sampling of it will be degenerate. Lawrence Murray PMCMC in marine biogeochemistry : 13 of 41
PZ model Simple phytoplankton-zooplankton model [Jones et al., 2010]. Lotka-Volterra with stochastic growth rate and quadratic mortality term. dP = α t P − cPZ dt dZ = ecPZ − m l Z − m q Z 2 dt α t ∼ N ( µ, σ ) . P is observed with noise. Simulated data used here. Lawrence Murray PMCMC in marine biogeochemistry : 14 of 41
PZ model: state estimate 7 Prior 2.5 Posterior Observed 2 6 Truth 1.5 Z 5 1 0.5 4 0 20 40 60 80 100 P 3 2 1 2 0 α 1 -1 -2 0 0 20 40 60 80 100 0 20 40 60 80 100 t t Lawrence Murray PMCMC in marine biogeochemistry : 15 of 41
NPZD model The model features nine noise terms, ξ i for i = 1 , . . . , 9 , each coupled to a univariate autoregressive process B i . B i ( t + ∆ t ) = B i ( 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 : 16 of 41
NPZD model A multiplicative temperature correction Tc is applied to all rate processes, for which a Q 10 formulation for dependence on temperature, T , is used: Tc = Q ( T − T ref ) / 10 , 10 where T ref is a reference temperature, and Q 10 a prescribed constant. Lawrence Murray PMCMC in marine biogeochemistry : 17 of 41
NPZD model The zooplankton grazing rate ( gr ) is dependent on the phytoplankton concentration (zooplankton functional response): gr = Tc · I Z · A υ (1 + A υ ) , (1) where υ is a given power. Lawrence Murray PMCMC in marine biogeochemistry : 18 of 41
NPZD model The relative availability of phytoplankton, A , is A = Cl Z · P ; I Z I Z is the maximum zooplankton ingestion rate (mg P per mg Z per day); Cl Z is the maximum clearance rate (volume in m 3 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 : 19 of 41
NPZD model A quadratic formulation for zooplankton mortality is adopted after Steele [1976] and Steele and Henderson [1992]: m = Tc · m Q · Z, where the quadratic mortality rate m Q has units of d − 1 ( mg Z m − 3 ) − 1 . Lawrence Murray PMCMC in marine biogeochemistry : 20 of 41
NPZD model The detrital remineralization rate is dependent only on temperature: r = Tc · r D , where r D prescribes the remineralisation rate at a reference temperature. Lawrence Murray PMCMC in marine biogeochemistry : 21 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 · g max · h E · h N / ( h E + h N ) . Lawrence Murray PMCMC in marine biogeochemistry : 22 of 41
NPZD model The light-limitation factor is given by h E = 1 − exp( − α · λ max · E/g max ) , where α is the initial slope of the photosynthesis versus irradiance curve (mg C mg Chla − 1 mol photon − 1 m 2 ), and λ max is the maximum Chla : C ratio (mg Chla mg C − 1 ). Lawrence Murray PMCMC in marine biogeochemistry : 23 of 41
NPZD model E is the mean photosynthetic available radiation (PAR) in the mixed layer and is given by E = E 0 . (1 − exp( − Kz )) /Kz, where E 0 is the mean daily photosynthetically available radiation (PAR) just below the air-sea interface and Kz is given by Kz = ( K W + a Ch · Chla ) · MLD. (2) K W is attenuation due to the seawater and a Ch . Lawrence Murray PMCMC in marine biogeochemistry : 24 of 41
NPZD model The nutrient-limitation factor is given by N h N = ( g max · Tc/a N ) + N , where a N is the maximum specific affinity for nitrogen uptake (d − 1 mg N − 1 m 3 ). Lawrence Murray PMCMC in marine biogeochemistry : 25 of 41
NPZD model The phytoplankton N : C ratio, χ , predicted by the model is given by χ = χ min · h E + χ max · h N , h E + h N where χ min and χ max are the prescribed minimum and maximum N : C ratios (mg N mg C − 1 ). Lawrence Murray PMCMC in marine biogeochemistry : 26 of 41
NPZD model The equations governing interactions between the remaining state variables { P, Z, D, N } are: dP κ = g · P − gr · Z + MLD · ( BCP − P ) dt dZ = E Z · gr · Z − m · Z dt dD D = (1 − E Z ) · f D · gr · Z + m · Z − r · D − S D · MLD + dt κ MLD · ( BCD − D ) dN = − g · P + (1 − E Z ) · (1 − f D ) · gr · Z + r · D + dt κ MLD · ( BCN − N ) . Lawrence Murray PMCMC in marine biogeochemistry : 27 of 41
State estimation (observed) 350 300 DIN ( µ g N l −1 ) 250 200 150 100 50 0 1971 1972 1973 1974 1975 4 ln(Chla ( µ g Chla l −1 )) 2 0 −2 −4 −6 −8 1971 1972 1973 1974 1975 Prior−95% Posterior−95% Observations Prior−median Posterior−median Lawrence Murray PMCMC in marine biogeochemistry : 28 of 41
Recommend
More recommend
