Discussion: Structured Forms for Model Discrepancy¶
Description and background¶
The basic ideas of model discrepancy are presented in the variant thread on linking models to reality (ThreadVariantModelDiscrepancy). There, the key equations are \(y=f(x^+) + d\) and \(z=y+\epsilon\), where \(y\) represents system values, \(f\) is the simultaor, \(x^+\) is the best input and \(z\) are observations on \(y\) with measurement error \(\epsilon\) (DiscObservations). We expand that notation here to allow explicit indexing of these quantities by other inputs; for example, space-time co-ordinates. Possible situations where \(d\) has a highly structured form due to physical structures present in the model are discussed.
Discussion¶
In accordance with the standard toolkit notation, we denote the model simulator by \(f\) and its inputs by \(x\). We extend the notation used in ThreadVariantModelDiscrepancy to include generalised ‘location’ inputs \(u\); for example, \(u\) might index space-time, such as latitude, longitude and date, and \(y(u)\) might denote the real average temperature of the earth’s surface at that latitude and longitude over a particular day. We write \(f(u, x)\) and \(y(u)\) to denote the simulator output for inputs \(x\) at location \(u\) and the system value at location \(u\), respectively. The fundamental relationship linking model to reality given in ThreadVariantModelDiscrepancy now becomes
where \(d(u)\) is the model discrepancy between the real system value \(y(u)\) at \(u\) and the simulator value \(f(u, x^+)\) for the best input \(x^+\) at \(u\): for detailed discussions about \(x^+\) and \(d(u)\) see DiscBestInput and DiscWhyModelDiscrepancy. As there may be replicate observations of the real system at a location \(u\), we introduce notation to cover this case. Denote by \(z_r(u)\) the \(r\)-th replicate observation of the system value \(y(u)\) at location \(u\). As in ThreadVariantModelDiscrepancy and DiscObservations, we assume that the \(z_r(u)\) and the real system value \(y(u)\) are additively related through the observation equation
where \(\epsilon_r(u)\) is the measurement error of the \(r\)-th replicate observation on \(y(u)\). We write \(\epsilon(u)\) when there is no replication. Combining these equations, we obtain the following relationship between system observations and simulator output
Statistical assumptions¶
- The components of the collection \(f(u, \cdot )\), \(x^+\), \(d(u)\), \(\epsilon_r(u)\) are assumed to be independent random quantities for each \(u\) in the collection \(U\) of \(u\)-values considered, or just uncorrelated in the Bayes linear case (ThreadCoreBL).
- The distribution of \(f(u, \cdot )\), \(x^+\), \(d(u)\), \(\epsilon_r(u)\) over all finite collections of \(u\)-values in \(U\) needs to be specified in any particular application.
- The model discrepancy term \(d(\cdot)\) is regarded as a random process over the collection \(U\). Either the distribution of \(d(\cdot)\) is specified over all finite subsets of \(U\), such as a Gaussian process for a full Bayes analysis (ThreadCoreGP), or just the expectation and covariance structure of \(d(\cdot)\) are specified over all finite subsets of \(U\) as is required for a Bayes linear analysis (ThreadCoreBL). A separable covariance structure for \(d(\cdot)\), similar to that defined in DefSeparable, might be a convenient simple choice when, for example, \(u\) indexes space-time.
- As the system is observed at only a finite number of locations \(u_1, \ldots, u_k\) (say), unlike \(d(\cdot)\), we do not regard \(\epsilon_r(\cdot)\) as a random process over \(U\). In most applications, the measurement errors \(\epsilon_{r_i}(u_i)\) are assumed to be independent and identically distributed random quantities with zero mean and variance \(\Sigma_\epsilon\). However, measurement errors may be correlated across the locations \(u_1, \ldots, u_k\), which might be the case, for example, when \(u\) indexes time. On the other hand, the replicate observations within each of locations \(u_1, \ldots, u_k\) are almost always assumed to be independent and identically distributed; see DiscObservations.
Additional comments and Examples¶
Kennedy, M. C. and O’Hagan, A. (2001) generalise the combined equation to
where \(\rho\) is a parameter to be specified or estimated.
Rainfall runoff¶
Iorgulescu, I., Beven, K. J., and Musy, A. (2005) consider a rainfall runoff model which simulates consecutive hourly measurements of water discharge \(D\) and Calcium \(C\) and Silica \(S\) concentrations for a particular catchment area. Here, \(u=t\) is hour \(t\) and \(y(t)= (D(t),C(t),S(t))\) is the amount of water discharged into streams and the Calcium and Silica concentrations at hour \(t\). There were 839 hourly values of \(t\). While the authors did not consider model discrepancy explicitly, a simple choice would be to specify the covariance \(\textrm{Cov}[d(t_k),d(t_{\ell})]\) between \(d(t_k)\) and \(d(t_{\ell})\) to be of the form
A more realistic simple choice might be the covariance structure derived from the autoregressive scheme \(d(t_{k+1})=\rho d(t_k) + \eta_k\), where \(\eta_k\) is a white noise process with variance \(\sigma^2\); the parameters \(\rho\) and \(\sigma^2\) are to be specified or estimated. Such a covariance structure would be particularly appropriate when forecasting future runoff.
The input vector \(x\) has eighteen components. An informal assessment of model discrepancy for this runoff model is given in Goldstein, M., Seheult, A. and Vernon, I. (2010): see also DiscInformalAssessMD.
Galaxy formation¶
Goldstein, M. and Vernon, I. (2009) consider the galaxy formation model ‘Galform’ which simulates two outputs, the \(b_j\) and \(K\) band luminosity functions. The \(b_j\) band gives numbers of young galaxies \(s\) per unit volume of different luminosities, while the \(K\) band describes the number of old galaxies \(l\). The authors consider 11 representative outputs, 6 from the \(b_j\) band and 5 from the \(K\) band. Here, \(u=(A,\lambda)\) is age \(A\) and luminosity \(\lambda\) and \(y(A,\lambda)\) is count of age \(A\) galaxies of luminosity \(\lambda\) per unit volume of space. The authors carried out a careful elicitation process with the cosmologists for \(\log y\) and specified a covariance \(\textrm{Cov}[d(A_i,\lambda_l),d(A_j,\lambda_k)]\) between \(d(A_i, \lambda_l)\) and \(d(A_j, \lambda_k)\) of the form
for specified values of the overall variance \(a\), the correlation within bands \(b\) and the correlation between bands \(c\). The input vector \(x\) has eight components.
Hydrocarbon reservoir¶
Craig, P. S., Goldstein, M., Rougier, J. C., and Seheult, A. H. (2001) consider a hydrocarbon reservoir model which simulates bottom-hole pressures of different wells through time. Here, \(u=(w,t)\) is the pair well \(w\) and time \(t\) and \(y(w,t)\) is the bottom-hole pressure for well \(w\) at time \(t\). There were 34 combinations of \(w\) and \(t\) considered. The authors specify a non-separable covariance \(\textrm{Cov}[d(w_i,t_k),d(w_j,t_{\ell})]\) between \(d(w_i,t_k)\) and \(d(w_j,t_{\ell})\) of the form
where \(I_P\) denotes the indicator function of the proposition \(P\). The input vector \(x\) has four active components: see DefActiveInput.
Spot welding¶
Higdon, D., Kennedy, M., Cavendish, J. C., Cafeo, J. A., and Ryne, R. D. (2004) consider a model for spot welding which simulates spot weld nugget diameters for different combinations of load, and current applied to two metal sheets, and gauge is the thickness of the two sheets. Here, \(u=(l,c,g)\) is the triple load \(l\), current \(c\) and gauge \(g\) and \(y(l, c, g)\) is the weld diameter when load \(l\) and current \(c\) are applied to sheets of gauge \(g\) at the \(12=2 \times 3 \times 2\) combinations. Moreover, there is system replication of 10 observations for each of the 12 system combinations. The authors specify a Gaussian process for the model discrepancy \(d\) over \((l, c, g)\) with a separable covariance structure. The input vector \(x\) has one component.
References¶
Craig, P. S., Goldstein, M., Rougier, J. C., and Seheult, A. H. (2001), “Bayesian forecasting for complex systems using computer simulators”, Journal of the American Statistical Association, 96, 717-729.
Goldstein, M. and Vernon, I. (2009), “Bayes linear analysis of imprecision in computer models, with application to understanding the Universe”, in 6th International Symposium on Imprecise Probability: Theories and Applications.
Goldstein, M., Seheult, A. and Vernon, I. (2010), “Assessing Model Adequacy”, MUCM Technical Report 10/04.
Higdon, D., Kennedy, M., Cavendish, J. C., Cafeo, J. A., and Ryne, R. D. (2004), “Combining field data and computer simulations for calibration and prediction”, SIAM Journal on Scientific Computing, 26, 448-466.
Iorgulescu, I., Beven, K. J., and Musy, A. (2005), “Data-based modelling of runoff and chemical tracer concentrations in the Haute-Mentue research catchment (Switzerland)”, Hydrological Processes, 19, 2557-2573.
Kennedy, M. C. and O’Hagan, A. (2001), “Bayesian calibration of computer models”, Journal of the Royal Statistical Society, Series B, 63, 425-464.