Computational Approaches for Parameter Estimation in Climate Models - PowerPoint PPT Presentation

Computational Approaches for Parameter Estimation in Climate Models Gabriel Huerta Department of Mathematics and Statistics University of New Mexico Joint with A. Villagran (UNM), C. Jackson and M.K. Sen (UT-Austin) 1

Summary Points • Parametric uncertainties in climate modeling. • Estimation of multidimensional probability distributions. • Climate model, proxy or surrogate to a full ACGM. • Simulated Annealing based method (MVFSA). • Adaptive Metropolis as an alternative. • Comparisons of posterior probability distributions.

Comments on Climate Models • IPCC Third Assessment Report (TAR, 2001): Global temperatures are likely to increase by 1.1 to 6.4 ◦ C 1990- 2100. • ”more comprehensive and systematic system of model analysis and diagnosis, and a Monte Carlo approach to model uncertainties associated with parameterizations” . • Recent progress with models of reduced complexity (For- est et al., 2000, 2001, 2002). • Perturbed physics ensembles with a general circulation model (Allen, 1999; Murphy et al, 2004; Stainforth et al., 2005; Collins et al., 2006).

Range of Model Hierarchies • General Circulation Models: Most demanding. – 16 processors/24 hours to simulate 10 years of cli- mate. • Models of Reduced Complexities: One or more spatial dimensions are eliminated. • Surrogate or Emulator Models: – Mimics equilibrium space-time response of an At- mospheric GCM. – To test sampling strategies for parametric uncertain- ties.

Goals of Present Study • Estimate probability distributions for parameters in cli- mate models. • Non-standard methods for state-of-art in the climate lit- erature. • Improve on the calibration of the climate model. • ”appreciate the ability to detect and attribute the effects of forcings on paleoclimate observations”. • Bayesian approach via posterior distributions.

Surrogate Climate Model • Jackson and Broccoli (2003): ”Changes in Earth’s orbital geometry over the past 165kyears”. • Obliquity Φ ∈ (22 ◦ , 25 ◦ ) , eccentricity, e ∈ (0 , 0 . 05) and longitude of perihelion, λ ∈ (0 ◦ , 360 ◦ ) . • d obs ( i, j, k ) is the observed surface temperature at lati- tude i , longitude j and season k . • Data simulated using Φ = 22 . 62 , e = 0 . 044 , λ = 75 . 93 . • d obs ( i, j, k ) = g ijk ( m ) + η ijk , where m = (Φ , e, λ ) . • g is the ”surface air temperature” to a given change in parameters.

• The surrogate model is g ijk ( m ) = A o Φ ′ + eA p cos ( φ p − λ ) + R ijk • Φ = Φ o + A o Φ ′ , Φ o is the mean obliquity. Φ ′ is the deviation of obliquity from its 165,000 year • mean. • φ p is the phase of the response to precessional forcing. • R ijk is a term that is added to represent the effects of internal variability. • A o and A p are the sensitivity of temperature to changes in obliquity and precession respectively.

Cost or Misfit Function • Measure of the deviation from the observed data and the model. • For the climate model considered, I J K E ( m ) = 1 � � � B − 1 ijk ( d obs ( i, j, k ) − g ijk ( m )) 2 2 N i =1 j =1 k =1 • B ijk is the variance of the observations at each grid point. • Other functions under study.

• The likelihood function is L ( d obs | m, S ) ∝ exp {− SE ( m ) } • S weights the significance of model-data differences. • We fixed S = 1 , S = 47 and plot a profile likelihood for each parameter. • Simulation values ( Φ = 22 . 625 , e = 0 . 043954 , λ = 75 . 93 ). • 20,000 point grid evaluation.

Likelihood functions for each orbital parameter. S=1 S=47 22 23 24 25 22 23 24 25 Obliquity 0 100 200 300 0 100 200 300 Longitude of Perihelion 0 0.01 0.02 0.03 0.04 0.05 0 0.01 0.02 0.03 0.04 0.05 Eccentricity

Multiple Very Fast Simulated Annealing • Ingber (1989). Given a current selection m ( k ) i , m ( k +1) = m ( k ) + y i ( m max − m min ) i i i i • y i ∼ Cauchy distribution. • Uses a cooling schedule T k = T o exp ( − α ( k − 1) 1 /d ) • One accepts or rejects m ( k +1) according to a Metropolis i scheme. • Sen and Stoffa (1996), multiple repetitions. • Balance between estimating a PPD and finding the global minimum.

Start at m and a given temperature T 0 Evaluate E( m ) 0 Optimization method Multiple VFSA Metropolis Is the number of perturbations NO ≤ moves/temperature ? Draw a number 0 ≤ y ≤ 1 YES Lower temperature, T from a white distribution Draw a number 0 ≤ y (T) ≤ 1 from a Cauchy distribution which is temperature, T, dependent Compute new model components m = m + y( m - m ) new 0 max min i i i i Evaluate E( m ) new YES NO Is D E=E( m ) - E( m ) < 0? Accept m with a new new 0 m = m 0 new probability exp(- D E/T) YES Has m changed in previous ntarget perturbations? STOP 0

Adaptive Methods: Single Component AM m i,k − 1 = ( m (0) i , ..., m ( k − 1) • ) are the sampled values for i the i-th parameter. • Variance equation V ( k ) = s d V ( m i,k − 1 ) + s d ǫ where i k − 1 1 ( m ( r ) � − m i ) 2 V ( m i,k − 1 ) = i k − 1 r =0 • Updating m ( k ) at iteration k , i – z i ∼ N ( m ( k − 1) , V ( k ) ) i i

– Accept z i with probability � 1 , π ( m ( k ) 1 , ..., m ( k ) i − 1 , z i , m ( k − 1) i +1 , ..., m ( k − 1) � ) d min π ( m ( k ) 1 , ..., m ( k ) i − 1 , m ( k − 1) , ..., m ( k − 1) ) i d • We produce samples π ( m, S | d obs ) . • ”Flat” priors on orbital forcing parameters. S ∼ Gamma distribution ( α 0 = 552 . 25 , β 0 = 11 . 75 ). • π ( m | S, d obs ) is sampled through SC-AM. π ( S | m, d obs ) is a Gamma α ∗ = α 0 and β ∗ = β 0 + • E ( m ( k ) ) .

Adaptive Methods: FAM. (Haario et al. 2001) • All parameters are sampled at once. • z ∼ q t ( ·| m (0) , ..., m ( t − 1) ) , • z is accepted with probability, � � π ( z ) α ( m ( t − 1) , z ) = min 1 , , π ( m ( t − 1) ) q t ( · ) is a multivariate Gaussian with mean m ( t − 1) and • covariance matrix C t . • Recursively, C t = s d C t − 1 + s d ǫI d , where s d > 0 and ǫ > 0

Adaptive Methods: DRAM (Haario and Mira 2006) • Combines Delayed Rejection (DR) and Adaptive Metropolis. • In DR, propose a second state move if first move was rejected. • The process can be iterated for a fixed or random num- ber. • DR combines different proposals. • May use for initial period (burn-in). • The DR method uses rejected values without losing prop- erties.

Bivariate scatter plots of orbital forcing parameters. FAM SCAM DRAM MVFSA METSA 300 300 300 300 300 Longitude of Perihelion 200 200 200 200 200 100 100 100 100 100 0 0 0 0 0 22 23 24 25 22 23 24 25 22 23 24 25 22 23 24 25 22 23 24 25 Obliquity Obliquity Obliquity Obliquity Obliquity 0.05 0.05 0.05 0.05 0.05 0.04 0.04 0.04 0.04 0.04 Eccentricity 0.03 0.03 0.03 0.03 0.03 0.02 0.02 0.02 0.02 0.02 0.01 0.01 0.01 0.01 0.01 0 0 0 0 0 22 23 24 25 22 23 24 25 22 23 24 25 22 23 24 25 22 23 24 25 Obliquity Obliquity Obliquity Obliquity Obliquity 0.05 0.05 0.05 0.05 0.05 0.04 0.04 0.04 0.04 0.04 Eccentricity 0.03 0.03 0.03 0.03 0.03 0.02 0.02 0.02 0.02 0.02 0.01 0.01 0.01 0.01 0.01 0 0 0 0 0 0 100 200 300 0 100 200 300 0 100 200 300 0 100 200 300 0 100 200 300 Longitude of Longitude of Longitude of Longitude of Longitude of Perihelion Perihelion Perihelion Perihelion Perihelion

Root Mean Square (RMS) Probability Error • For every parameter θ , : RMS i ( θ ) = || Prob ( θ ) − Prob ( θ ) π || i • Prob ( θ ) is a vector of ”bin probabilities” estimated from i output at iteration i . • Prob ( θ ) probability vector under target . i • RMS goes to zero as i goes to infinity. • Prob ( θ ) from a baseline method. i

Comparison for Obliquity Parameter. Quantile Estimation − 97.5% 1.8 25 FAM SCAM 1.6 DRAM MVFSA 24.5 METSA 1.4 1.2 24 Obliquity 1 23.5 0.8 0.6 23 FAM 0.4 SCAM DRAM 22.5 MVFSA 0.2 METSA 0 22 22 22.5 23 23.5 24 24.5 25 0 10000 20000 30000 40000 50000 Obliquity Iterations Obliquity 25 25 FAM SCAM 24.5 DRAM 20 MVFSA METSA 24 15 Obliquity RMS 23.5 10 23 5 22.5 0 22 0 10000 20000 30000 40000 50000 FAM SCAM DRAM MVFSA METSA Iterations

Comparison for Longitude of Perihelion Parameter. Quantile Estimation − 97.5% 0.06 350 FAM SCAM DRAM 0.05 300 MVFSA METSA Longitude of Perihelion 250 0.04 FAM SCAM 200 DRAM 0.03 MVFSA METSA 150 0.02 100 0.01 50 0 0 0 50 100 150 200 250 300 350 0 10000 20000 30000 40000 50000 Longitude of Perihelion Iterations Longitude of Perihelion 0.35 350 FAM SCAM 0.3 300 DRAM MVFSA 0.25 METSA Longitude of Perihelion 250 0.2 200 RMS 0.15 150 0.1 100 0.05 50 0 0 0 10000 20000 30000 40000 50000 FAM SCAM DRAM MVFSA METSA Iterations