Model Uncertainty: Mathematical and Statistical

By E. Bruce Pitman

Mathematical models developed for computational simulation of complex, real-world processes are crucial ingredients in virtually every field of science, engineering, medicine, and business. The remarkable growth in computing power and the matching gains in algorithmic speed and accuracy have increased the applicability and reliability of simulations, making them faster and enabling their employment in previously intractable problems. Statistical models play a crucial role in analyzing real-world processes in virtually every field. The recent explosion in the storage and availability of data—and the computing power required to handle it—has increased the centrality of statistical analysis. A convenient way of categorizing models is to say that statistical models are useful for processes with abundant available data, and math models are valuable in data-poor environments.

It is rarely—perhaps never—possible to construct a mathematical or statistical model of a process with assurance (from its construction) that the model accurately represents all of the process’s details. The crucial final step when utilizing a model for prediction or comprehension of a real-world process is to understand the uncertainties inherent in that model’s use. Mathematicians and engineers call this process uncertainty quantification, which has become a major part of applied mathematics. Statisticians refer to it as model uncertainty (MU) — one of the most prominent fields of statistics.

In 2018-2019, the Statistical and Applied Mathematical Sciences Institute (SAMSI) hosted a year-long program called “Model Uncertainty: Mathematical and Statistical” (MUMS), with several workshops occurring throughout the year. SAMSI’s central research activity involves working groups: collections of faculty visitors, postdoctoral fellows, and students with a shared interest in a particular topic or problem. The MUMS program generated both working groups with a more theoretical leaning and working groups of a more applied nature, and sponsored activities that fostered exchanges between the theoretical and applied. Programming activities in statistics, modeling, and computational mathematics illustrate MUMS’ success and highlight some important problems in MU.

Mathematical modelers commonly formulate a large computer simulation of a physical or biological process, composed of a discretization of a system of ordinary differential equations (ODEs) and/or partial differential equations (PDEs). Meaningful data with which to judge a model requires fusion of results from computer simulations and experimental measurements. This approach gives rise to several questions. To begin, can one calibrate the model and determine good values for the important parameters in the system that are consistent with the data? Answering this question is predicated on determining which parameters are important — itself a difficult task. Next comes calibration, a procedure that typically involves a Markov chain Monte Carlo (MCMC) method to calculate the probability distribution for these parameters. MCMC demands many evaluations of the governing differential equations, meaning that constructing a surrogate for the governing system is often helpful.

In statistics, Gaussian Stochastic Process (GaSP) emulators comprise the workhorse surrogate model. GaSP emulators are a nonparametric regression of data that are equipped with an estimate of the likely error of that regression. But when it comes to model behavior at a large number of space/time points, the typical GaSP construction is computationally intractable, thus necessitating new approaches. Efficient construction of a GaSP emulator for high-dimensional inputs and outputs—and the associated problem of finding reduced-order models for the system of interest—touched many of the working groups and represented a recurring theme throughout the year.

The preceding discussion makes a fundamental assumption: that the computed system of equations adequately represents the physical or biological processes being studied. A moment’s reflection tells us that no mathematical model can capture every relevant process across all length and time scales. Thus, a somewhat more precise question is: Given a specified quantity of interest, does the mathematical model being computed satisfactorily represent the operative processes? That is, is the model system adequate for purpose?

Bob Moser (University of Texas at Austin) presented new ideas on this problem of model discrepancy to one of the working groups, which stimulated months-long discussions. Perhaps a simple example can illustrate the important questions at the heart of the challenge: the motion of a damped linear spring (see Figure 1). Moser asked participants to imagine a damped spring with a mass at its end, the wrinkle being that the constitutive relation for the damping is unknown. After receiving data consisting of the position of the mass at several times for mass values 1 and 2, one is asked to predict the velocity of the spring mass at a specified time for mass value 5. This problem presents two subtle issues. First, the data gives the position value at several times, but researchers are required to predict a latent variable: the velocity. Second, data is provided for values of the mass equal to 1 and 2, but the question seeks to extrapolate to mass 5.

Figure 1. Data for the damped spring experiment. Red and green dots show the position of the mass at times \(t=1, ... 8\) for mass 1 and mass 2 respectively. The blue curve, which is not known to the analyst, indicates the position for mass \(=5\); the challenge is to predict the velocity at a specified time. Figure courtesy of E. Bruce Pitman and Pierre Barbillon.

Let us examine this problem a bit more deeply. Given the data, we model the mass’s location at the \(i\)th time as 

\[ z_i = \zeta(t_i) + \epsilon_i,\]

where \(\zeta\) is the true value of the spring extension for the correct input parameters and \(\epsilon\) represents observation errors, assumed to be independent and normally distributed. A computational model of the spring system encodes conservation of mass and momentum with initial and boundary conditions. However, the constitutive relation between stress and stress rate—and strain and strain rate—is at best an approximation based on prior experiments. We believe that the conservation principles hold universally but know the constitutive law holds only in a limited range of conditions. It is here that model discrepancy comes into play [2-4]. We postulate an error due to model inadequacy (owing to the unknown damping relation):

\[\zeta(t) =\eta(t, \theta) + \delta(t). \tag1 \]

Here, \(\theta\) is a set of calibration parameters that are not exactly known. Thus, our observation model is

\[z_i =\eta(t_i, \theta) + \delta(t_i)  + \epsilon_i. \tag2\]

That is, the observation is a sum of the “true” value of the system, an error due to model inadequacy, and measurement error. The classic Kennedy and O’Hagan framework [4] models \(\delta\) as a GaSP, assigns a prior to \(\theta\), and determines posterior distributions conditioned on the observed data.

Moser described ideas published in [5], wherein a discrepancy operator is postulated and embedded in the dynamical system. An earlier example of this idea contains the damped spring-mass scenario [6]. Because the actual damping is unknown, including model inadequacy attempts to rectify this constitutive ignorance. The embedded discrepancy operator contains several parameters that must also be calibrated, adding to the burden of analysis.

An alternative approach to representing model error provides moment-based likelihoods through a generalized polynomial chaos expansion [7]. Assuming that the form of model discrepancy continues to hold, the chaos expansion can be valuable when attempting to make predictions for parameters outside of range testing.

A different approach to studying dynamic discrepancy considers the suspect parameters as varying in time and adds a dynamic discrepancy: \(\delta (x, \chi(t))\). One may assume a Gaussian process prior, which turns the governing ODEs or PDEs into a system of stochastic differential equations [1]. Dave Mebane (West Virginia University) presented this idea at the 2019 Materials and Data Science Hackathon, sponsored by SAMSI and the Data-enabled Science and Engineering of Atomic Structure group at North Carolina State University.

Imagine a system of two equations that take the form

\[\frac{\partial x}{\partial t} = v;\]

\[ \frac{\partial v}{\partial t} = - c(\chi)  v - k x  + \delta({ x},  \chi(t), \beta),\]

where \(\beta\) are hyperparameters of the Gaussian process. If more than one parameter is suspect, expanding them all can explode the number of terms to be determined. Therefore, a low-order Karhunen-Loève or polynomial chaos expansion is often specified to reduce the work.

Pierre Barbillon (AgroParisTech) illustrated yet another approach at the MUMS Transition Workshop. In our joint work, he extends the governing system by postulating a stochastic relaxation. The governing equations take the form of a generalized mean-reverting process with a random forcing. The spring-mass system becomes

\[\frac{\partial x}{\partial t} = v; \:\:\:\:\:\:\:\:\: m \frac{\partial v}{\partial t} = w\]

\[ \textrm{and}\]

\[\frac{\partial w}{\partial t} = -\frac{1}{\tau} (w - (-c v  - kx)) + \sigma^2 \frac{\partial B}{\partial t}.\]

Here we abuse notation in writing \(\partial B/\partial t,\), where \(B\) is Brownian motion. Notice that in the zero-relaxation limit \(\tau \rightarrow 0\), the auxiliary variable \(w \rightarrow -(cv+kx)\) formally, thus recovering the original ODEs. Simulations allow this system to wander away from the (inadequate) constant coefficient linear spring; calibration \(\tau\) and \(\sigma^2\) suggest the magnitude of this discrepancy.

This example of model inadequacy suggests important issues of model and variable selection (the focus of another working group). It highlights the need for a close collaboration between domain scientists, statisticians, and computational and applied mathematicians to evaluate MU in a way that allows quantification of the uncertainty in predictions, which is the investigation’s ultimate goal.

References
[1] Bhat, K.S., Mebane, D.S., Storlie, C.B., & Mahapatra, P. (2017). Upscaling uncertainty with dynamic discrepancy for a multi-scale carbon capture system. J. Amer. Stat. Assoc., 112, 1453-1467.
[2] Brynjarsdottir, J., & O’Hagan, A. (2014). Learning about physical parameters: The importance of model discrepancy. Inv. Prob., 30(11).
[3] O’Hagan, A. (2013). Bayesian inference with mis-specied models: Inference about what? J. Stat. Plan. Infer., 143, 1643-1648.
[4] Kennedy, M.C., & O’Hagan, A. (2001). Bayesian calibration of computer models. J. Roy. Soc. Stats. B, 63, 425-464.
[5] Morrison, R., Oliver, T., & Moser, R.D. (2016). Representing model inadequacy: A stochastic operator approach. SIAM J. Uncert. Quant., 6, 457-496.
[6] Oliver, T., Terejanu, G., Simmons, C.S., & Moser, R.D. (2015). Validating predictions of unobserved quantities. Comp. Meth. Appl. Mech. Eng., 283, 1310-1335.
[7] Sargsyan, K., Huan, X., & Njam, H.N. (2018). Embedded model error representation for Bayesian model calibration. Int. J. Uncert. Quant.

E. Bruce Pitman is a professor in the Department of Materials Design and Innovation at the University at Buffalo. He chaired the organizing committee for the 2018-2019 “Model Uncertainty: Mathematical and Statistical” program at the Statistical and Applied Mathematical Sciences Institute.