Search

Many dynamical models, such as numerical weather prediction and climate models, contain so called closure parameters. These parameters usually appear in physical parameterizations of sub-grid scale processes, and they act as "tuning handles" of the models. Currently, the values of these parameters are specified mostly manually, but the increasing complexity of the models calls for more algorithmic ways to perform the tuning. Traditionally, parameters of dynamical systems are estimated by directly comparing the model simulations to observed data using, for instance, a least squares approach. However, if the models are chaotic, the classical approach can be ineffective, since small errors in the initial conditions can lead to large, unpredictable deviations from the observations. In this paper, we study numerical methods available for estimating closure parameters in chaotic models. We discuss three techniques: off-line likelihood calculations using filtering methods, the state augmentation method, and the approach that utilizes summary statistics from long model simulations. The properties of the methods are studied using a modified version of the Lorenz 95 system, where the effect of fast variables are described using a simple parameterization.

The Kalman filter (KF) and Extended Kalman filter (EKF) are well-known tools for assimilating data and model predictions. The filters require storage and multiplication of $n\times n$ and $n\times m$ matrices and inversion of $m\times m$ matrices, where $n$ is the dimension of the state space and $m$ is dimension of the observation space. Therefore, implementation of KF or EKF becomes impractical when dimensions increase. The earlier works provide optimization-based approximative low-memory approaches that enable filtering in high dimensions. However, these versions ignore numerical issues that deteriorate performance of the approximations: accumulating errors may cause the covariance approximations to lose non-negative definiteness, and approximative inversion of large close-to-singular covariances gets tedious. Here we introduce a formulation that avoids these problems. We employ L-BFGS formula to get low-memory representations of the large matrices that appear in EKF, but inject a stabilizing correction to ensure that the resulting approximative representations remain non-negative definite. The correction applies to any symmetric covariance approximation, and can be seen as a generalization of the Joseph covariance update. &nbsp We prove that the stabilizing correction enhances convergence rate of the covariance approximations. Moreover, we generalize the idea by the means of Newton-Schultz matrix inversion formulae, which allows to employ them and their generalizations as stabilizing corrections.

In a typical inverse problem, a spatially distributed parameter in a physical model is estimated from measurements of model output. Since measurements are stochastic in nature, so is any parameter estimate. Moreover, in the Bayesian setting, the choice of regularization corresponds to the definition of the prior probability density function, which in turn is an uncertainty model for the unknown parameters. For both of these reasons, significant uncertainties exist in the solution of an inverse problem. Thus to fully understand the solution, quantifying these uncertainties is important. When the physical model is linear and the error model and prior are Gaussian, the posterior density function is Gaussian with a known mean and covariance matrix. However, the electrical impedance tomography inverse problem is nonlinear, and hence no closed form expression exists for the posterior density. The typical approach for such problems is to sample from the posterior and then use the samples to compute statistics (s...

The extended Kalman filter (EKF) is a popular state estimation method for nonlinear dynamical models. The model error covariance matrix is often seen as a tuning parameter in EKF, which is often simply postulated by the user. In this paper, we study the filter likelihood technique for estimating the parameters of the model error covariance matrix. The approach is based on computing the likelihood of the covariance matrix parameters using the filtering output. We show that (a) the importance of the model error covariance matrix calibration depends on the quality of the observations, and that (b) the estimation approach yields a well-tuned EKF in terms of the accuracy of the state estimates and model predictions. For our numerical experiments, we use the two-layer quasi-geostrophic model that is often used as a benchmark model for numerical weather prediction.

High-dimensional inverse problems present a challenge for Markov chain Monte Carlo (MCMC)-type sampling schemes. Typically, they rely on finding an efficient proposal distribution, which can be dif...

In the ensemble Kalman filter (EnKF), uncertainty in the state of a dynamical model is represented as samples of the state vector. The samples are propagated forward using the evolution model, and the forecast (prior) mean and covariance matrix are estimated from the ensemble. Data assimilation is carried out by using these estimates in the Kalman filter formulas. The prior is given in the subspace spanned by the propagated ensemble, the size of which is typically much smaller than the dimension of the state space. The rank-deficiency of these covariance matrices is problematic, and, for instance, unrealistic correlations often appear between spatially distant points, and different localization or covariance tapering methods are needed to make the approach feasible in practice. In this paper, we present a novel way to implement ensemble Kalman filtering using optimization-based sampling, in which the forecast error covariance has full rank and the need for localization is diminished. The method is based on the randomize then optimize (RTO) technique, where a sample from a Gaussian distribution is computed by perturbing the data and the prior, and solving a quadratic optimization problem. We test our method in two benchmark problems: the 40dimensional Lorenz’96 model and the 1600-dimensional two-layer quasi-Geostrophic model. Results show that the performance of the method is significantly better than that of the standard EnKF, especially with small ensemble sizes when the rank-deficiency problems in EnKF are emphasized.

Algorithmic numerical weather prediction (NWP) skill optimization has been tested using the Integrated Forecasting System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF). We report the results of initial experimentation using importance sampling based on model parameter estimation methodology targeted for ensemble prediction systems, called the ensemble prediction and parameter estimation system (EPPES). The same methodology was earlier proven to be a viable concept in low-order ordinary differential equation systems, and in large-scale atmospheric general circulation models (ECHAM5). Here we show that prediction skill optimization is possible even in the context of a system that is (i) of very high dimensionality, and (ii) carefully tuned to very high skill. We concentrate on four closure parameters related to the parameterizations of sub-grid scale physical processes of convection and formation of convective precipitation. We launch standard ensembles of medium-range predictions such that each member uses different values of the four parameters, and make sequential statistical inferences about the parameter values. Our target criterion is the squared forecast error of the 500 hPa geopotential height at day three and day ten. The EPPES methodology is able to converge towards closure parameter values that optimize the target criterion. Therefore, we conclude that estimation and cost function-based tuning of low-dimensional static model parameters is possible despite the very high dimensional state space, as well as the presence of stochastic noise due to initial state and physical tendency perturbations. The remaining question before EPPES can be considered as a generally applicable tool in model development is the correct formulation of the target criterion. The one used here is, in our view, very selective. Considering the multi-faceted question of improving forecast model performance, a more general target criterion should be developed. This is a topic of ongoing research.

The reaction kinetics of styrene and indene with cyclopentadiene in the presence of tert-butylcatechol as a polymerization inhibitor was studied in a batch reactor in the production of potential monomers for amorphous polymers: 5-phenyl-2-norbornene and indenylnorbornene. The experimental results can be explained with a reaction mechanism using four separate reactions: an equilibrium reaction of dicyclopentadiene and cyclopentadiene, a Diels–Alder reaction of cyclopentadiene with a dienophile to produce norbornene, disproportionation of norbornene, and a reaction between norbornene and dicyclopentadiene. Some polymerization could be noticed when an excess of dienophile was used. The reaction kinetic parameters of this mechanism were fitted and their reliabilities tested using Markov chain Monte Carlo (MCMC) methods. The suggested reaction mechanisms explained the experimental results well, but analysis with MCMC methods revealed clear uncertainties in the reliabilities of the estimated parameters. Here, M...

The ensemble Kalman filter (EnKF) is a technique for dynamic state estimation. EnKF approximates the standard extended Kalman filter (EKF) by creating an ensemble of model states whose mean and empirical covariance are then used within the EKF formulas. The technique has a number of advantages for large-scale, nonlinear problems. First, large-scale covariance matrices required within EKF are replaced by low rank and low storage approximations, making implementation of EnKF more efficient. Moreover, for a nonlinear state space model, implementation of EKF requires the associated tangent linear and adjoint codes, while implementation of EnKF does not. However, for EnKF to be effective, the choice of the ensemble members is extremely important. In this paper, we show how to use the conjugate gradient (CG) method, and the recently introduced CG sampler, to create the ensemble members at each filtering step. This requires the use of a variational formulation of EKF. The effectiveness of the method is demonstrated on both a large-scale linear, and a small-scale, non-linear, chaotic problem. In our examples, the CG-EnKF performs better than the standard EnKF, especially when the ensemble size is small.

Prior distributions for Bayesian inference that rely on the $l_1$-norm of the parameters are of considerable interest, in part because they promote parameter fields with less regularity than Gaussian priors (e.g., discontinuities and blockiness). These $l_1$-type priors include the total variation (TV) prior and the Besov $B^s_{1,1}$ space prior, and in general yield non-Gaussian posterior distributions. Sampling from these posteriors is challenging, particularly in the inverse problem setting where the parameter space is high-dimensional and the forward problem may be nonlinear. This paper extends the randomize-then-optimize (RTO) method, an optimization-based sampling algorithm developed for Bayesian inverse problems with Gaussian priors, to inverse problems with $l_1$-type priors. We use a variable transformation to convert an $l_1$-type prior to a standard Gaussian prior, such that the posterior distribution of the transformed parameters is amenable to Metropolized sampling via RTO. We demonstrate this approach on several deconvolution problems and an elliptic PDE inverse problem, using TV or Besov $B^s_{1,1}$ space priors. Our results show that the transformed RTO algorithm characterizes the correct posterior distribution and can be more efficient than other sampling algorithms. The variable transformation can also be extended to other non-Gaussian priors., Comment: Preprint 24 pages, 13 figures. v1 submitted to SIAM Journal on Scientific Computing on July 5, 2016

We apply a recently developed method, the Ensemble Prediction and Parameter Estimation System (EPPES), to demonstrate how numerical weather prediction (NWP) model closure parameters can be optimized. As proof of concept, we tune the medium-range forecast skill of the ECMWF model HAMburg version (ECHAM5) atmospheric general circulation model using an ensemble prediction system (EPS) emulator. Initial state uncertainty is represented in the EPS emulator by applying the initial state perturbations generated at the European Centre for Medium-range Weather Forecasts (ECMWF). Model uncertainty is represented in the emulator via parameter variations at the initial time. We vary four closure parameters related to parametrizations of subgrid-scale physical processes of clouds and precipitation. With this set-up, we generate ensembles of 10-day global forecasts with the ECHAM5 model at T42L31 resolution twice a day over a period of three months. The cost function in the optimization is formulated in terms of standard forecast skill scores, verified against the ECMWF operational analyses. A summarizing conclusion of the experiments is that the EPPES method is able to find ECHAM5 model closure parameter values that correspond to smaller values of the cost function. The forecast skill score improvements verify positively in dependent and independent samples. The main reason is the reduced temperature bias in the tropical lower troposphere. Moreover, the optimization improved the top-of-atmosphere radiation flux climatology of the ECHAM5 model, as verified against the Clouds and the Earth's Radiant Energy System (CERES) radiation data over a 6-year period, while the simulated tropical cloud cover was reduced, thereby increasing a negative bias as verified against the International Satellite Cloud Climatology Project (ISCCP) data.

In model-based process optimization one uses a mathematical model to optimize a certain criterion, for example the product yield of a chemical process. Models often contain parameters that have to be estimated from data. Typically, a point estimate (e.g. the least squares estimate) is used to fix the model for the optimization stage. However, parameter estimates are uncertain due to incomplete and noisy data. In this article, it is shown how parameter uncertainty can be taken into account in process optimization. To quantify the uncertainty, Markov Chain Monte Carlo (MCMC) sampling, an emerging standard approach in Bayesian estimation, is used. In the Bayesian approach, the solution to the parameter estimation problem is given as a distribution, and the optimization criteria are functions of that distribution. The formulation and implementation of the optimization is studied, and numerical examples are used to show that parameter uncertainty can have a large effect in optimization results.

Summary 1. World-wide extinctions of amphibians are at the forefront of the biodiversity crisis, with climate change figuring prominently as a potential driver of continued amphibian decline. As in other taxa, changes in both the mean and variability of climate conditions may affect amphibian populations in complex, unpredictable ways. In western North America, climate models predict a reduced duration and extent of mountain snowpack and increased variability in precipitation, which may have consequences for amphibians inhabiting montane ecosystems. 2. We used Bayesian capture–recapture methods to estimate survival and transition probabilities in a high-elevation population of the Columbia spotted frog (Rana luteiventris) over 10 years and related these rates to interannual variation in peak snowpack. Then, we forecasted frog population growth and viability under a range of scenarios with varying levels of change in mean and variance in snowpack. 3. Over a range of future scenarios, changes in mean snowpack had a greater effect on viability than changes in the variance of snowpack, with forecasts largely predicting an increase in population viability. Population models based on snowpack during our study period predicted a declining population. 4. Although mean conditions were more important for viability than variance, for a given mean snowpack depth, increases in variability could change a population from increasing to decreasing. Therefore, the influence of changing climate variability on populations should be accounted for in predictive models. The Bayesian modelling framework allows for the explicit characterization of uncertainty in parameter estimates and ecological forecasts, and thus provides a natural approach for examining relative contributions of mean and variability in climatic variables to population dynamics. 5. Longevity and heterogeneous habitat may contribute to the potential for this amphibian species to be resilient to increased climatic variation, and shorter-lived species inhabiting homogeneous ecosystems may be more susceptible to increased variability in climate conditions.

Classical optimal design criteria suffer from two major flaws when applied to nonlinear problems. First, they are based on linearizing the model around a point estimate of the unknown parameter and therefore depend on the uncertain value of that parameter. Second, classical design methods are unavailable in ill-posed estimation situations, where previous data lack the information needed to properly construct the design criteria. Bayesian optimal design can, in principle, solve these problems. However, Bayesian design methods are not widely applied, mainly due to the fact that standard implementations for efficient and robust routine use are not available. In this article, we point out a concrete recipe for implementing Bayesian optimal design, based on the concept of simulation-based design introduced by Muller, Sanso, and De Iorio (2004). We develop further a predictive variance criterion and introduce an importance weighting mechanism for efficient computation of the variances. The simulation-based appr...

SUMMARY The Kalman filter is a technique for estimating a time-varying state given a dynamical model for and indirect measurements of the state. It is used, for example, on the control problems associated with a variety of navigation systems. Even in the case of nonlinear state and/or measurement models, standard implementations require only linear algebra. However, for sufficiently large-scale problems, such as arise in weather forecasting and oceanography, the matrix inversion and storage requirements of the Kalman filter are prohibitive, and hence, approximations must be made. In this paper, we describe how the conjugate gradient iteration can be used within the Kalman filter for quadratic minimization, as well as for obtaining low-rank approximations of the covariance and inverse-covariance matrices required for its implementation. The approach requires that we exploit the connection between the conjugate gradient and Lanczos iterations. Copyright © 2011 John Wiley & Sons, Ltd.

In the companion paper (Jarvinen et al., 2011. Q. J. R. Meteorol. Soc. DOI: 10.1002/qj.923) we suggested a concept to estimate numerical weather prediction model closure parameters on-line, coupled with an operational ensemble prediction system. In this paper, such a method is developed and demonstrated in low-order numerical tests. The method utilizes the massive amount of operational model predictions to make statistical inference about the underlying probability distributions of the closure parameters. This is achieved with practically no additional computations. The only required change to the ensemble prediction system is to allow perturbation of the model closure parameters for different ensemble members. Otherwise, the method is straightforward. The parametric uncertainty is presented using a hierarchical statistical model. Proposed parameter values are resampled based on their respective likelihood function values, as evaluated against verifying observations. Update formulas are used to feed back the parametric information to the evolving proposal distribution. The method for ensemble prediction and parameter estimation system (EPPES) is demonstrated using a stochastic version of the Lorenz-95 model. The numerical tests show that the EPPES method is capable of detecting unknown and wrongly specified parameter values, and lead to optimal forecast skill in an independent test run. Potential show-stoppers are discussed. Our current research is directed towards a realistic ensemble prediction system. Copyright © 2011 Royal Meteorological Society

This paper discusses algorithmic estimation of numerical weather prediction (NWP) model closure parameters which appear in physical parametrizations of sub-grid scale processes. Currently, these parameters are specified manually. Our hypothesis is that the NWP model forecast skill is sensitive to the specified parameter values. Therefore some parameter value combinations score better than others, but it is very demanding to manually specify the optimal combination. First, we approximate the number of test forecasts that would be needed to estimate the posterior probability densities of a limited number of closure parameters, in a spirit of the Bayesian parameter estimation approach. Second, we discuss various options to accommodate the related huge computational task; we end up proposing that operational ensemble prediction might be a potential host. Third, we analyse the ensemble prediction system characteristics from the parameter estimation point of view. We rule out various standard parameter estimation methods; instead we suggest an estimation concept that can be on-line coupled with existing ensemble prediction systems to utilize the operational ensemble simulations. Thus practically no additional computations are introduced. Such a method, the ensemble prediction and parameter estimation system (EPPES), is developed here and demonstrated in low-order numerical tests in the companion paper (Laine M, Solonen A, Haario H, Jarvinen H. 2011. Q. J. R. Meteorol. Soc. DOI: 10.1002/qj.922). Copyright © 2011 Royal Meteorological Society

Climate models contain closure parameters to which the model climate is sensitive. These parameters appear in physical parameterization schemes where some unresolved variables are expressed by predefined parameters rather than being explicitly modeled. Currently, best expert knowledge is used to define the optimal closure parameter values, based on observations, process studies, large eddy simulations, etc. Here, parameter estimation, based on the adaptive Markov chain Monte Carlo (MCMC) method, is applied for estimation of joint posterior probability density of a small number (n=4) of closure parameters appearing in the ECHAM5 climate model. The parameters considered are related to clouds and precipitation and they are sampled by an adaptive random walk process of the MCMC. The parameter probability densities are estimated simultaneously for all parameters, subject to an objective function. Five alternative formulations of the objective function are tested, all related to the net radiative flux at the top of the atmosphere. Conclusions of the closure parameter estimation tests with a low-resolution ECHAM5 climate model indicate that (i) adaptive MCMC is a viable option for parameter estimation in large-scale computational models, and (ii) choice of the objective function is crucial for the identifiability of the parameter distributions.

Synthesis of the oxygenate fuel additive TAME (tertiary amyl methyl ether) is based on the reaction between reactive isopentene compounds (2-methyl-1-butene and 2-methyl-2-butene) and methanol. Skeletal isomerization of n-pentenes is an important reaction route in production of the raw material for TAME synthesis. In this work the kinetics of skeletal isomerization of 1-pentene was studied using HZSM-22 as a catalyst, and the effect of temperature (523–598 K) and WHSV (weight hourly space velocity) on product composition were investigated. The observed main reaction products were cis- and trans-2-pentenes, 2-methyl-2-butene, 2-methyl-1-butene and 3-methyl-1-butene. Minor products found in the output stream were C2, C3, C4, C6, C7, C8 and C10 – alkenes and – alkanes. Only traces of methane and C9 compounds were detected. Two different reaction mechanisms have been proposed in literature for skeletal isomerization of n-pentenes and their side reactions (dimerization and fragmentation); the monomolecular and the bimolecular skeletal isomerization mechanisms. In the current work a pseudohomogeneous reaction model for both mechanisms is fitted to experimental data. The fitting results show that both reaction mechanisms explain the skeletal isomerization of 1-pentene to isopentenes well. However, the monomolecular reaction mechanism explained the formation of dimerization and fragmentation products slightly better. The reliability of the estimated Arrhenius parameters was analyzed using Bayesian inference and Markov chain Monte Carlo methods (MCMC) which are rarely used in chemical engineering.

The reaction kinetics of p-nitrophenol ozonation in a bubble column at pH 2 of water is studied. The reaction kinetic parameters, rate coefficients and stoichiometric coefficients, as well as the mass transfer coefficient are estimated by nonlinear estimation with the developed novel reaction kinetic model. The concentration of unknown intermediate compounds is presented as a residual chemical oxygen demand (COD) calculated from the measured COD and the theoretical COD for the known compounds. The decomposition rate of p-nitrophenol is almost the same on the pathway producing hydroquinone and on the pathway producing 4-nitrocathecol. Comparison of the rate coefficients with the case of initial pH 5 indicates that the p-nitrophenol degradation producing 4-nitrocathecol is more selective towards molecular ozone than the reaction producing hydroquinone. The identifiability and reliability of the estimated parameters, as well as the uncertainties involved in the model predictions are analyzed with the Markov chain Monte Carlo (MCMC) method. MCMC methods are Bayesian statistical methods that can be used for analyzing the distribution of parameters in nonlinear models in order to extend and support the information given by traditional regression analysis.

The aim of this work was to study the reaction kinetics of neopentyl glycol esterification with different carboxylic acids (propionic acid, isobutyric acid and 2-ethylhexanoic acid). The experiments were carried out in a batch reactor as reactive distillation removing continuously one reaction product (water) by purging nitrogen through the reactor. Commercial heterogeneous gellular (Dowex 50WX2) and macroreticular (Amberlyst 15) type resins and homogeneous para-toluene sulfonic acid were compared as catalysts. The gellular type resin was dried before loading into the reactor. When performing esterification of neopentyl glycol with different acids the observed formation of monoester was clearly much faster than the consecutive esterification of monoester to diester. A disproportionation reaction was also found to take place. The reliability of the estimated Arrhenius parameters was analyzed using Bayesian inference and Markov chain Monte Carlo methods (MCMC). MCMC analysis was also used to study uncertainties in the model predictions and in some additional parameters, such as reaction rate constants calculated from the estimated parameters. This cannot be done by classical regression analysis.

The reaction kinetics of neopentyl glycol esterification with propionic acid is studied in a batch reactor at different temperatures. The reaction is catalyzed by a heterogeneous resin (DOWEX 50 WX 2). The parameters of the kinetic model are estimated by non-linear regression by minimizing the difference between the measured and predicted concentrations. The reliability of the estimated parameters is analyzed using the Markov chain Monte Carlo (MCMC) method. The MCMC method is an advanced statistical approach, which is well suited for the analysis of complex multiparametric non-linear models. The potential of the MCMC method is further demonstrated by optimizing the product yield in CSTR using the estimated parameters from the batch experiments.

In this work, the reaction kinetics of p-nitrophenol ozonation in a bubble column at low pH of water was studied. The reaction kinetic parameters, rate coefficients, and stoichiometric coefficients were estimated by using nonlinear optimization. The reaction rate equations were written by applying a part of the reaction scheme of Yu and Yu (Ozone Sci. Eng. 2001, 23, 303−312). The concentration of the unknown intermediate compounds is presented as a residual COD (chemical oxygen demand) calculated from the measured COD and the theoretical COD for the known species. The decomposition rate of p-nitrophenol on the pathway producing hydroquinone was found to be about two times faster than the decomposition rate on the pathway producing 4-nitrocatechol. The ozone self-decomposition by hydroperoxide ions seems to be the plausible cause for the relatively high ozone consumption at low pH. The identifiability and reliability of the estimated parameters were analyzed with the Markov chain Monte Carlo (MCMC) method....

In this work, we consider the Bayesian optimization (BO) approach for tuning parameters of complex chaotic systems. Such problems arise, for instance, in tuning the sub-grid scale parameterizations in weather and climate models. For such problems, the tuning procedure is generally based on a performance metric which measures how well the tuned model fits the data. This tuning is often a computationally expensive task. We show that BO, as a tool for finding the extrema of computationally expensive objective functions, is suitable for such tuning tasks. In the experiments, we consider tuning parameters of two systems: a simplified atmospheric model and a low-dimensional chaotic system. We show that BO is able to tune parameters of both the systems with a low number of objective function evaluations and without the need of any gradient information.

In the Bayesian approach to inverse problems, data are often informative, relative to the prior, only on a low-dimensional subspace of the parameter space. Significant computational savings can be achieved by using this subspace to characterize and approximate the posterior distribution of the parameters. We first investigate approximation of the posterior covariance matrix as a low-rank update of the prior covariance matrix. We prove optimality of a particular update, based on the leading eigendirections of the matrix pencil defined by the Hessian of the negative log-likelihood and the prior precision, for a broad class of loss functions. This class includes the F\"{o}rstner metric for symmetric positive definite matrices, as well as the Kullback-Leibler divergence and the Hellinger distance between the associated distributions. We also propose two fast approximations of the posterior mean and prove their optimality with respect to a weighted Bayes risk under squared-error loss. These approximations are deployed in an offline-online manner, where a more costly but data-independent offline calculation is followed by fast online evaluations. As a result, these approximations are particularly useful when repeated posterior mean evaluations are required for multiple data sets. We demonstrate our theoretical results with several numerical examples, including high-dimensional X-ray tomography and an inverse heat conduction problem. In both of these examples, the intrinsic low-dimensional structure of the inference problem can be exploited while producing results that are essentially indistinguishable from solutions computed in the full space.

The intrinsic dimensionality of an inverse problem is affected by prior information, the accuracy and number of observations, and the smoothing properties of the forward operator. From a Bayesian perspective, changes from the prior to the posterior may, in many problems, be confined to a relatively low-dimensional subspace of the parameter space. We present a dimension reduction approach that defines and identifies such a subspace, called the 'likelihood-informed subspace' (LIS), by characterizing the relative influences of the prior and the likelihood over the support of the posterior distribution. This identification enables new and more efficient computational methods for Bayesian inference with nonlinear forward models and Gaussian priors. In particular, we approximate the posterior distribution as the product of a lower-dimensional posterior defined on the LIS and the prior distribution marginalized onto the complementary subspace. Markov chain Monte Carlo sampling can then proceed in lower dimensions, with significant gains in computational efficiency. We also introduce a Rao−Blackwellization strategy that de-randomizes Monte Carlo estimates of posterior expectations for additional variance reduction. We demonstrate the efficiency of our methods using two numerical examples: inference of permeability in a groundwater system governed by an elliptic PDE, and an atmospheric remote sensing problem based on Global Ozone Monitoring System (GOMOS) observations., United States. Dept. of Energy. Office of Advanced Scientific Computing Research (Grant DE-SC0003908), United States. Dept. of Energy. Office of Advanced Scientific Computing Research (Grant DE-SC0009297)

Parameterisation schemes of subgrid-scale physical processes in atmospheric models contain so-called closure parameters. Their precise values are not generally known; thus, they are subject to fine-tuning for achieving optimal model performance. In this article, we show that there is a dilemma concerning the optimal parameter values: an identical prediction model formulation can have two different optimal closure parameter value settings depending on the level of approximations made in the data assimilation component of the prediction system. This result tends to indicate that the prediction model re-tuning in large-scale systems is not only needed when the prediction model undergoes a major change, but also when the data assimilation component is updated. Moreover, we advocate an accurate albeit expensive method based on so-called filter likelihood for the closure parameter estimation that is applicable in fine-tuning of both prediction model and data assimilation system parameters. In this article, we use a modified Lorenz-95 system as a prediction model and extended Kalman filter and ensemble adjustment Kalman filter for data assimilation. With this setup, we can compute the filter likelihood for the chosen parameters using the output of the two versions of the Kalman filter and apply a Markov chain Monte Carlo algorithm to explore the parameter posterior distributions. Keywords: model tuning, Markov chain Monte Carlo, likelihood, filter formulation (Published: 3 May 2013) Citation: Tellus A 2013, 65 , 20147, http://dx.doi.org/10.3402/tellusa.v65i0.20147

A priori dimension reduction is a widely adopted technique for reducing the computational complexity of stationary inverse problems. In this setting, the solution of an inverse problem is parameterized by a low-dimensional basis that is often obtained from the truncated Karhunen–Loève expansion of the prior distribution. For high-dimensional inverse problems equipped with smoothing priors, this technique can lead to drastic reductions in parameter dimension and significant computational savings. In this paper, we extend the concept of a priori dimension reduction to non-stationary inverse problems, in which the goal is to sequentially infer the state of a dynamical system. Our approach proceeds in an offline–online fashion. We first identify a low-dimensional subspace in the state space before solving the inverse problem (the offline phase), using either the method of 'snapshots' or regularized covariance estimation. Then this subspace is used to reduce the computational complexity of various filtering algorithms—including the Kalman filter, extended Kalman filter, and ensemble Kalman filter—within a novel subspace-constrained Bayesian prediction-and-update procedure (the online phase). We demonstrate the performance of our new dimension reduction approach on various numerical examples. In some test cases, our approach reduces the dimensionality of the original problem by orders of magnitude and yields up to two orders of magnitude in computational savings., Academy of Finland (Projects 284715 and 267442), United States. Dept. of Energy. Office of Advanced Scientific Computing Research (Grant DE-SC0003908)

The emergence of Markov chain Monte Carlo (MCMC) methods has opened a way for Bayesian analysis of complex models. Running MCMC samplers typically requires thousands of model evaluations, which can exceed available computer resources when this evaluation is computationally intensive. We will discuss two generally applicable techniques to improve the efficiency of MCMC. First, we consider a parallel version of the adaptive MCMC algorithm of Haario et al. (2001), implementing the idea of inter-chain adaptation introduced by Craiu et al. (2009). Second, we present an early rejection (ER) approach, where model simulation is stopped as soon as one can conclude that the proposed parameter value will be rejected by the MCMC algorithm. ¶ This work is motivated by practical needs in estimating parameters of climate and Earth system models. These computationally intensive models involve non-linear expressions of the geophysical and biogeochemical processes of the Earth system. Modeling of these processes, especially those operating in scales smaller than the model grid, involves a number of specified parameters, or ‘tunables’. MCMC methods are applicable for estimation of these parameters, but they are computationally very demanding. Efficient MCMC variants are thus needed to obtain reliable results in reasonable time. Here we evaluate the computational gains attainable through parallel adaptive MCMC and Early Rejection using both simple examples and a realistic climate model.

Studying different theoretical properties of epidemiological models has been widely addressed, while numerical studies and especially the calibration of models, which are often complicated and loaded with a high number of unknown parameters, against measured data have received less attention. In this paper, we describe how a combination of simulated data and Markov Chain Monte Carlo (MCMC) methods can be used to study the identifiability of model parameters with different type of measurements. Three known models are used as case studies to illustrate the importance of parameter identifiability: a basic SIR model, an influenza model with vaccination and treatment and a HIV–Malaria co-infection model. The analysis reveals that calibration of complex models commonly studied in mathematical epidemiology, such as the HIV–Malaria co-dynamics model, can be difficult or impossible, even if the system would be fully observed. The presented approach provides a tool for design and optimization of real-life field campaigns of collecting data, as well as for model selection.