capbtabboxtable[][]

# Orthogonal calibration via posterior projections with applications to the
Schwarzschild model

Antik Chakraborty<EMAIL_ADDRESS>Department of Statistics, Purdue
University  Jonelle L. Walsh<EMAIL_ADDRESS>George P. and Cynthia
Woods Mitchell Institute for Fundamental Physics and Astronomy, Department of
Physics and Astronomy, Texas A&M University  Louis Strigari<EMAIL_ADDRESS>George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and
Astronomy, Department of Physics and Astronomy, Texas A&M University  Bani K.
Mallick<EMAIL_ADDRESS>Department of Statistics, Texas A&M University
Anirban Bhattacharya<EMAIL_ADDRESS>Department of Statistics, Texas
A&M University

###### Abstract

The orbital superposition method originally developed by Schwarzschild (1979)
is used to study the dynamics of growth of a black hole and its host galaxy,
and has uncovered new relationships between the galaxy’s global
characteristics. Scientists are specifically interested in finding optimal
parameter choices for this model that best match physical measurements along
with quantifying the uncertainty of such procedures. This renders a
statistical calibration problem with multivariate outcomes. In this article,
we develop a Bayesian method for calibration with _multivariate outcomes_
using orthogonal bias functions thus ensuring parameter identifiability. Our
approach is based on projecting the posterior to an appropriate space which
allows the user to choose any nonparametric prior on the bias function(s)
instead of having to model it (them) with Gaussian processes. We develop a
functional projection approach using the theory of Hilbert spaces. A finite-
dimensional analogue of the projection problem is also considered. We
illustrate the proposed approach using a BART prior and apply it to calibrate
the Schwarzschild model illustrating how a multivariate approach may resolve
discrepancies resulting from a univariate calibration.

Keywords: Bayesian; Computer model; Emulator; Scientific modeling; Regression
trees.

##  1 Introduction

Scientific endeavors typically aim to describe a natural physical phenomenon.
In modern applications, these descriptions increasingly involve complex
mathematical equations typically solved using expensive computer codes. For
example, the MIT2D climate model
(http://web.mit.edu/globalchange/www/climate.html) simulates probability
distributions of ocean, surface and upper atmospheric temperatures. In our
motivating application, Schwarzschild’s orbital integral method
(Schwarzschild, 1979) is used to simulate mass distributions of galaxies.
These orbit-based models are the main method for measuring the masses of
supermassive black holes at the centers of nearby galaxies (Kormendy and Ho,
2013), and have led to the establishment of empirical relationships between
the masses of black holes and large-scale galaxy properties, which
surprisingly suggest that black holes and their host galaxies grow and evolve
together over time (McConnell and Ma, 2013; Kormendy and Ho, 2013; Saglia et
al., 2016). Besides determining black hole masses, Schwarzschild modeling is a
powerful tool for measuring a galaxy’s mass-to-light ratio, dark matter halo
properties, stellar orbital distribution, viewing orientation, and intrinsic
three-dimensional shape, allowing for the further study of galaxy assembly
histories (Mehrgan et al., 2019). These scientific models involve parameters
which carry physical meaning. Inferring about the parameters while addressing
the uncertainty of the data observed from the the physical process was first
formulated as a statistical problem in Sacks et al. (1989), and was further
studied from a Bayesian point of view in Kennedy and O’Hagan (2001), who
popularized the term computer code calibration. Since then statistical
calibration methods have been used in many scientific disciplines; e.g.
climatology (Forest et al., 2008; Salter et al., 2018), cell biology (Xie and
Xu, 2021), mechanical engineering (Gattiker et al., 2006) etc.

Parallel to the widening applications of statistical calibration methods,
substantial effort has also been dedicated to the development of more robust
calibration methods within the statistics community, especially after the
seminal work of Kennedy and O’Hagan (2001). In their work, the authors
proposed a framework that simultaneously models the data observed from a
physical process and data generated from computer code implementations of
mathematical models emulating the physical process. The fundamental
contribution of their work was to include a discrepancy/bias function in their
model, which the authors interpreted as a unknown function that captures the
inability of the mathematical model to explain variations in the observed
data. A Bayesian nonparametric regression approach based on Gaussian processes
(Rasmussen and Williams, 2005) was then used to model the unknown bias
function as well as the computer code simulating the physical process. This
enabled fast emulation of the computer code simulators and produced rapid
estimates of input parameters of the computer code that best approximated the
observed data. Since then several authors have proposed statistical methods
for computer code calibration. Higdon et al. (2008) developed a basis vector
approach to calibration when the computer code output is very high-dimensional
using principal components; Bayarri et al. (2007) proposed a general framework
for Bayesian calibration where a modular Bayes approach was advocated to
overcome significant computational and inferential challenges associated with
a fully Bayesian analysis of calibration problems; Chakraborty et al. (2013)
modeled the computer code using multivariate adaptive regression splines
(Friedman, 1991; Denison et al., 1998) to avoid overfitting the computer
model; Pratola and Higdon (2016) proposed to model the bias and computer model
using Bayesian additive regression trees (BART; Chipman et al. (2010)) which
has gained immense popularity in the machine learning community recently. Tuo
and Wu (2015, 2016) studied the general calibration model of Kennedy and
O’Hagan (2001) from a theoretical point-of-view and showed that the
calibration parameters are not identifiable. As a remedy, Tuo and Wu (2015)
defined the calibration parameter as the point in the input parameter space
that minimizes the $L^{2}$ loss function between the physical process and the
computer code. Plumlee (2017) extended it to more general losses and then
proposed a modified covariance function for Gaussian process priors used for
the bias function; Xie and Xu (2021) treated the calibration parameter as a
functional of the bias function and used a projected Bayesian approach to
avoid identifiability issues.

In this article, we propose a generalized Bayesian approach for computer code
calibration with multivariate output which takes into account identifiability
issues pointed out by Tuo and Wu (2016). We start from the definition of the
optimal calibration parameter as the point in the input space of the computer
code that minimizes a certain loss function. Following Plumlee (2017), this
identifiability issue can be mitigated by imposing a orthogonality constraint
on the bias function with respect to the derivatives of the computer code. We
first generalize the restrictions on the bias function to the situation where
there are multiple potentially correlated outcomes. In the univariate outcome
case, Plumlee (2017) developed a modified Gaussian process prior on the bias
function that respects the orthogonality constraint. However, as indicated by
other authors (Pratola and Higdon, 2016) Gaussian process priors often suffer
from scalability issues to big data sets and high-dimensional input spaces.
The scalability issues are only exacerbated when one has to deal with
multivariate covariance functions. On the other hand, while nonparametric
priors such as the BART are highly scalable, introducing a priori constraints
such as the orthogonality constraint indicated above is not immediate. To
overcome these issues, we propose to posit a unconstrained prior on the bias
function and project the obtained posterior on the relevant space. For
inference based on the projected posterior we develop two algorithms for
posterior sampling. The first one is essentially a Hilbert projection on the
space of square-integrable functions and the second one borrows inspiration
from a projection interpretation of the multivariate Gaussian distribution.
Numerical experiments in Section 5 show superiority of the proposed method
both in terms of computational scalability and modeling flexibility. Indeed,
the method provides a more general framework for orthogonal calibration
whereby practitioners are free to choose a nonparametric prior of their own
choice; for example deep Gaussian process priors used recently in Marmin and
Filippone (2022). We also note that our approach is significantly different
from a related projection based method by Xie and Xu (2021) where the authors
consider a Gaussian process prior on the bias function. The calibration
parameter, defined as the minimizer of an appropriate loss function between
the observed data and the computer model is then treated as a functional of
the bias function. Hence, a prior is induced on the calibration parameter
through the Gaussian process prior on the bias function. Although this serves
the statistical issue of identifiability, properties of the induced prior is
less understood.

The rest of the paper is organized as follows. In Section 2, we briefly
describe our motivating example which is then followed by the development of
the orthogonal posterior projection approach in Section 3. We then compare the
efficiency of the proposed method with other state-of-the-art methods in terms
of estimating the calibration parameter and computational scalability. We end
with an application on the computer code data from Schwazschild model. An R
package for implementing the proposed method is available to download from
https://github.com/antik015/OCal.

##  2 Schwarzschild’s model

The study of the mass distribution within galaxies is central in the quest of
understanding of black holes, stellar components, dark matter, and the growth
of galaxies. Schwarzschild (1979) presented a orbit superposition method for
constructing a self-consistent mass model of galaxies. It consists of
“integrating a representative set of orbits in a static triaxial gravitational
potential, and finding weights for these orbits such that their superposition
reproduces the assumed mass distribution” (Quenneville et al., 2021). Several
modifications to the initial method have since been proposed among which van
den Bosch et al. (2008)’s triaxial orbit superposition has become very
popular. In standard applications of the method, the model is compared with
kinematic and photometric data to determine best-fit parameters such as black
hole mass $(\theta_{1})$, stellar mass-to-light ratio $(\theta_{2})$, and
fraction of dark matter halo $(\theta_{2})$.

(a) Left: $v_{S}$; Right: $v_{F}$

(b) Left: $\tau_{S}$; Right: $\tau_{F}$

(c) Left: $h_{3,S}$; Right: $h_{3,F}$

(d) Left: $h_{4,S}$; Right: $h_{4,F}$

Figure 1: Orbital superposition output versus observed data on four moments of
the velocity distribution for
$\theta=(10.22,10,1)^{\mathrm{\scriptscriptstyle{T}}}$. The superposition
method is implemented by a code originally developed by van den Bosch et al.
(2008).

In it’s current version (van den Bosch et al., 2008), for a given input
parameter combination
$t=(t_{1},t_{2},t_{3})^{\mathrm{\scriptscriptstyle{T}}}$, the code outputs the
first four moments of the line-of-sight velocity distribution
$f_{S}(x;t)=(v_{S}(x;t),\tau_{S}(x;t),h_{3,S}(x;t),h_{4,S}(x;t))^{\mathrm{\scriptscriptstyle{T}}}$
for points $x$ in a 2-dimensional spatial grid $\mathcal{X}$. The points $x$
are typically chosen to match the locations of the physical photometric data,
providing further information about the mass distribution. Specifically,
physical data
$y_{F}(x_{i})=(v_{F}(x_{i}),\tau_{F}(x_{i}),h_{3,F}(x_{i}),h_{4,F}(x_{i}))^{\mathrm{\scriptscriptstyle{T}}}$
is available for $i=1,\ldots,n,\,n=105$ locations of $\mathcal{X}=[-1,1]^{2}$.
In Figure 1, color coded output from the code and the observed data are shown
when $t=(10.22,10,1)^{\mathrm{\scriptscriptstyle{T}}}$ at the $n$ locations.
This superposition technique is implemented by a computer code (FORTRAN) which
is extremely computationally intensive; e.g. for one input of $\theta$,
approximately 3 hours is needed to generate the output. The code output is
available for three different values of $t_{3}=1,2,3$, and for each value
$t_{3}$, the output is available over a two-dimensional grid of values for
$(t_{1},t_{2})$ where $t_{1}$ ranges from 10.22 to 10.31 and $t_{2}$ ranges
from 10 to 10.25. Overall, the computer code output is available for $N=476$
different values of $t$. Naturally, an exhaustive search over the input
parameter space that compares the code output to the observed data to find the
best-fit parameter combination is prohibitive. However, from limited
experiments performed separately on each of the 4 outputs of the model,
scientists have seen that for the point $\tilde{\theta}=(9.92,9.29,1.27)$ the
model outputs best approximated the observed data up to slight variations for
each of the outcome. It is, however, unclear how a joint analysis of the model
outputs and the physical data can be carried out that adequately addresses the
uncertainty in the observed data, and whether such an analysis would reveal
regions in the parameter space that perform relatively well. This is important
for designing future experiments where scientists plan to include other
parameters that could potentially lead to better predictions of the physical
process on previously unobserved points - a task typically suited for a
multivariate analysis when outcomes are correlated. Inspired by this
application, in Section 3, we first develop a general methodology for
orthogonal calibration of multivariate computer models where the key focus is
on parameter identifiability and modeling flexibility. The validity of the
proposed methodology is then investigated in a series of numerical experiments
in Section 5, followed by an application on the Schwarzschild’s model in
Section 6.

##  3 Method

###  3.1 Orthogonal calibration

We begin with a quick review of orthogonal calibration, and introduce notation
necessary to develop our method. Suppose we are interested in a real physical
system $y_{R}(\cdot)\in\mathbb{R}^{q},\,q\geq 1$ with controllable inputs
$x\in\mathcal{X}$ where $\mathcal{X}$ is some closed, bounded subset of
$\mathbb{R}^{d}$. We do not directly observe the system. Instead, we observe
stochastic field observations $y_{F}(\cdot)$ at $n$ input points called the
design $\mathcal{D}=\\{x_{1},\ldots,x_{n}\\}$, for which a natural stochastic
model is

$y_{F}(x_{i})=y_{R}(x_{i})+\epsilon(x_{i}),\,\epsilon(x_{i})\overset{ind.}{\sim}{\mathrm{N}}_{q}(0,\Sigma_{F}),\,i=1,\ldots,n.$
(1)

Suppose in addition to the field observations, we also have a simulator
$f:\mathcal{X}\times\Theta\to\mathbb{R}^{q}$ for the physical system
(depending on the controllable input $x$ as well as other parameters
$t\in\Theta\subset\mathbb{R}^{p}$) which typically comes in the form of a
computer code/model. The parameter $t$ represents our understanding of the
physical system, and we hope that for some unknown $\theta^{*}\in\Theta$, the
computer model closely approximates the physical process $y_{R}(\cdot)$.
Formally, suppose $L\\{y_{R}(\cdot),f(\cdot,\cdot)\\}$ is a loss function that
can distinguish between the best possible parameter $\theta^{*}$ and all other
parameter values $t$. If there exists a unique $t^{*}$ such that
$y_{R}(x)=f(x,t^{*})$ for all $x\in\mathcal{X}$ and the loss $L$ is strictly
proper following Gneiting and Raftery (2007), then $\theta^{*}=t^{*}$.
However, existence of such a $t^{*}$ is not always guaranteed for most
practical computer models and loss functions. In the absence of such a
parameter, the computer model is biased and $\theta^{*}$ is defined as the
parameter combination at which the loss $L\\{y_{R}(\cdot),f(\cdot,\cdot)\\}$
is minimized , i.e.
$\theta^{*}=\operatorname*{arg\,min}_{t\in\Theta}L\\{y_{R}(\cdot),f(\cdot,t)\\}$.
To account for such bias, Kennedy and O’Hagan (2001) proposed the following
model for the observed data on the system

$y_{F}(x_{i})=f(x_{i},\theta)+b_{\theta}(x_{i})+\epsilon(x_{i}),\,\epsilon(x_{i})\overset{ind.}{\sim}{\mathrm{N}}(0,\Sigma_{F}),\,i=1,\ldots,n,$
(2)

where $\theta$ represents the model parameter that targets the population
parameter $\theta^{*}$ defined above, and $b_{\theta}(x)$ is interpreted as
the discrepancy or bias between the physical process and the assumed computer
model, i.e. $b_{\theta}(x)=y_{R}(x)-f(x,\theta)$. Without added restrictions
in (2), $\theta$ is not identifiable in (2). We provide a sufficient condition
for identifiablity under the loss

$L\\{y_{R}(\cdot),f(\cdot,t)\\}=\sum_{k=1}^{q}\int_{\mathcal{X}}\\{y_{R,k}(x)-f_{k}(x,t)\\}^{2}dx,$
(3)

in the following proposition, extending the result of Plumlee (2017) to
multivariate outcomes.

######  Proposition 3.1.

Consider the loss $L\\{y_{R}(\cdot),f(\cdot,t)\\}$ defined above. Let
$g_{j,k}(x,t)=\frac{\partial}{\partial t_{j}}f_{k}(x,t_{j})$, the partial
derivative with respect to $t_{j}$ of the $k$-th component $f_{k}(x,t)$ of the
computer model $f(x,t)$, which is assumed to exist for all $j,k$. A sufficient
condition for $\frac{\partial}{\partial
t}L\\{y_{R}(\cdot),f(\cdot,t)\\}\rvert_{t=\theta^{*}}=0$ is

$\int_{\mathcal{X}}\sum_{k=1}^{q}g_{j,k}(x,\theta^{*})b_{\theta^{*},k}(x)dx=0,\quad\text{for
all }j=1,\ldots,p.$ (4)

###### Proof.

Write $y_{R,k}(x)=f_{k,\theta}(x)+b_{\theta,k}(x)$ to obtain that the $j$-th
element of

$\frac{\partial}{\partial
t}L\\{y_{R}(\cdot),f(\cdot,t)\\}\rvert_{t=\theta}=\\\
-2\int_{\mathcal{X}}\sum_{k=1}^{q}g_{j,k}(x)b_{\theta,k}(x)dx,$

and recall
$\theta^{*}=\operatorname*{arg\,min}_{t\in\Theta}L\\{y_{R}(\cdot),f(\cdot,t)\\}$
which completes the proof. ∎

In other words, if the bias functions satisfy a linear constraint with respect
to the gradients of the computer model at $\theta$, then $\theta$ is
identifiable.

###  3.2 Bayesian inference on model (2)

The influential paper Kennedy and O’Hagan (2001) originally proposed the model
(2) for $q=1$, and considered Gaussian process priors (Rasmussen and Williams,
2005) as a prior distribution over the unknown function $b_{\theta}(x)$
together with a Normal/Uniform prior over the parameter $\theta$. For now, we
assume that the functional form of the computer model $f(x,t)$ is explicitly
known for every $x\in\mathcal{X},t\in\Theta$ and defer the extension of the
proposed method to unavailable $f$ later. A Gaussian process prior over an
unknown function $h(\cdot)$ is typically elicited by a mean function
$m(\cdot)$ and a covariance function $C(\cdot,\cdot)$. Under this prior, for
any finite collection of points $\\{u_{1},\ldots,u_{N}\\}$, the prior
distribution over the vector
$\\{h(u_{1}),\ldots,h(u_{N})\\}^{\mathrm{\scriptscriptstyle{T}}}$ has a
multivariate Gaussian distribution with mean $\\{m(u_{1}),\ldots,m(u_{N})\\}$
and covariance matrix $K=(K_{ij})$ and $K_{ij}=C(u_{i},u_{j}),\,1\leq i,j\leq
N$. Let $\Pi(\theta)$ be the prior distribution over $\theta$ and suppose
$b_{\theta}$ is endowed with a Gaussian process prior with a zero mean
function and a covariance kernel $C(\cdot,\cdot)$. Then with $n$ observations
on $y_{F}(\cdot)$ on the design $\mathcal{D}$, the posterior distribution of
$(\theta,\mathbf{b}),\,\mathbf{b}=\\{b_{\theta}(x_{1}),\ldots,b_{\theta}(x_{n})\\}^{\mathrm{\scriptscriptstyle{T}}}$
can be obtained using Bayes theorem as

$\Pi(\theta,\mathbf{b}\mid
y^{(n)})\propto\ell(y^{(n)}\mid\theta,\mathbf{b})\Pi(\theta)\Pi(\mathbf{b}),$
(5)

where given $\theta$,
$y^{(n)}=\\{y_{F}(x_{1})-f(x_{1},\theta),\ldots,y_{F}(x_{n})-f(x_{n},\theta)\\}^{\mathrm{\scriptscriptstyle{T}}}$
and $\ell(y^{(n)}\mid\theta,\mathbf{b})$ is the likelihood corresponding to a
Gaussian distribution with mean $\mathbf{b}$ and covariance matrix
$\sigma^{2}_{F}\mathrm{I}_{n}$. The marginal posterior distribution of
$\theta$, $\Pi(\theta\mid y^{(n)})$ can be obtained by observing that
$y^{(n)}\mid\theta\sim{\mathrm{N}}(0,K+\sigma_{F}^{2}\mathrm{I}_{n})$.

Tuo and Wu (2016) proved that vanilla Gaussian process priors over the bias
$b_{\theta}(\cdot)$ may lead to inaccurate inference on $\theta^{*}$, defined
as the minimizer of the loss function discussed above, even when observations
are available for a large number of points $n$. This is mainly due to the fact
that a generic Gaussian process on the bias function is not guaranteed to
satisfy the orthogonality condition (4). In light of this, Plumlee (2017)
advocated a modified covariance kernel for $b_{\theta}(\cdot)$ such that
realizations of such a Gaussian process prior satisfy (4) almost surely, see
also Plumlee and Joseph (2018). This modified covariance kernel is

$C_{\theta}(x,x^{\prime})=C(x,x^{\prime})-h_{\theta}(x)^{\mathrm{\scriptscriptstyle{T}}}H_{\theta}^{-1}h_{\theta}(x^{\prime}),$
(6)

where the $j$-th element of the vector $h_{\theta}(x)\in\mathbb{R}^{p}$ is
given by $\int_{\mathcal{X}}g_{j}(x,\theta)C(x,u)du$ and the $(j,k)$-th
element of the matrix $H_{\theta_{0}}\in\mathbb{R}^{p\times p}$ is
$\int_{\mathcal{X}}\int_{\mathcal{X}}g_{j}(u,\theta)g_{k}(u^{\prime},\theta)C(u,u^{\prime})dudu^{\prime}$
with $x,x^{\prime}\in\mathcal{X}$ for $1\leq j,k\leq p$. This covariance
function is obtained by observing that under a Gaussian process prior on
$b_{\theta_{0}}$, the prior on
$(b_{\theta},\int_{\mathcal{X}}g_{1}(x,\theta)b_{\theta}(x)dx,\ldots,\int_{\mathcal{X}}g_{p}(x,\theta_{0})b_{\theta}(x)dx)$
is again a Gaussian process and one then simply finds the conditional
covariance of $b_{\theta}$ given that the last $p$ elements of the above
vector are 0, and thus the orthogonality condition in (4) is enforced. Of
course, apriori we do not have knowledge on $\theta$. However, assuming the
knowledge of the system $y_{R}(x)$ and the computer model $f(x,t)$ for all
$t\in\Theta$, a consistent estimator of $\theta$ can be used for implementing
orthogonal covariance function in practice (Plumlee, 2017).

###  3.3 Posterior projections

Although the modified covariance kernel seemingly resolves idenitifiability
issue of $\theta$, it hinges on the assumption that the user chooses to model
$b_{\theta}(x)$ with a Gaussian process. This leaves out a plethora of other
nonparametric priors for functions that have been widely popular in the
literature, for example Bayesian additive regression trees (BART) (Chipman et
al., 2010), Bayesian multivariate adaptive regression splines (Denison et al.,
1998), Bayesian neural networks (Neal, 2012) etc. Furthermore, these methods
have been successfully applied in computer code calibration problems (Pratola
and Higdon, 2016; Higdon et al., 2008; Chakraborty et al., 2013). Many of
these priors have been shown to have optimal theoretical properties along with
sufficient computational scalability (Ročková and van der Pas, 2020; Linero
and Yang, 2018). In this article, our aim is to develop a procedure where a
user is free to choose any nonparametric prior for $b_{\theta}(x)$ and the
procedure automatically takes care of orthogonality constraints discussed
above.

Let $L^{2}_{q}(\mathcal{X})$ be the space of $q$-dimensional square-integrable
functions on $\mathcal{X}$ equipped with the inner-product $\langle
f_{1},f_{2}\rangle=\sum_{k=1}^{q}\int f_{1,k}(x)f_{2,k}(x)dx$. Then
$L^{2}_{q}(\mathcal{X})$ is a tensor Hilbert space with norm
$\|\cdot\|_{L^{2}_{q}}$ induced by the inner product:
$\|f\|_{L^{2}_{q}}=\langle f,f\rangle=\sum_{k=1}^{q}\int f_{k}^{2}(x)dx$. We
assume that for every $t\in\Theta$,
$g_{j}(\cdot,t)=(g_{j,1}(\cdot,t),\ldots,g_{j,q}(\cdot,t))^{\mathrm{\scriptscriptstyle{T}}}\in
L^{2}_{q}(\mathcal{X})$ for all $j=1,\ldots,p$ and the corresponding bias
$b_{t}(\cdot)\in L^{2}_{q}(\mathcal{X})$. Define
$\mathcal{F}_{\tilde{\theta}}=\\{b_{\tilde{\theta}}:\int\sum_{k=1}^{q}g_{j,k}(x,\tilde{\theta})b_{\tilde{\theta},k}(x)dx=0,j=1,\ldots,p\\}$
for a fixed $\tilde{\theta}$. Here $\tilde{\theta}$ is our guess about the
parameter combination at which the computer model best approximates the system
given the observed data and the computer model. Such a data-dependent choice
of $\tilde{\theta}$ involves minimizing the corresponding empirical risk. Let
$\tilde{\theta}$ be defined as

$\tilde{\theta}=\operatorname*{arg\,min}_{t\in\Theta}\frac{1}{nq}\sum_{k=1}^{q}\sum_{i=1}^{n}\\{y_{F,k}(x_{i})-f_{k}(x_{i},t)\\}^{2}.$
(7)

Validity of this definition is provided in the following result where we first
establish $\theta^{*}$ as the minimizer of population risk and then leverage
on empirical risk minimization results to show that $\tilde{\theta}$ converges
to $\theta^{*}$ for large $n$.

######  Proposition 3.2.

Let $P=P_{y\mid x}P_{x}$ denote the data generating distribution where $P_{x}$
is the uniform measure on $\mathcal{X}$. Suppose that $\Theta$ is compact and
that the $k$-th computer model is Lipshitz in $t\in\Theta$ uniformly over
$x\in\mathcal{X}$, that is there exists $L_{k}>0$ such that
$|f_{k}(x,t_{1})-f_{k}(x,t_{2})|\leq L_{k}\|t_{1}-t_{2}\|_{2}^{2}$ for
$k=1,\ldots,q$. Then $\tilde{\theta}\overset{P}{\to}\theta^{*}$.

A proof is provided in the appendix. Our definition of $\tilde{\theta}$ is
different from the definition used in (Plumlee, 2017, Proposition 1) where the
author assumes the knowledge of $y_{R}(\cdot)$ along with the computer model
$f(\cdot,t)$ for all $t$. We argue that explicit functional form of
$y_{R}(\cdot)$ may not be available in many practical scenarios, and hence
define $\tilde{\theta}$ as detailed above. The compactness of $\Theta$ is a
technical assumption which nevertheless is often practical in many calibration
problems.

Having defined the point at which we want the bias function $b(x)$ to satisfy
the orthogonality condition (4), we now proceed to define the projection
posterior. Instead of defining a prior distribution on
$b_{\tilde{\theta}}\in\mathcal{F}_{\tilde{\theta}}$, our strategy is to start
with a generic prior and then project the posterior distribution to
$\mathcal{F}_{\tilde{\theta}}$. Standard Hilbert space theory implies that
$\mathcal{F}_{\tilde{\theta}}$ is a non-empty, closed and convex subset of
$L^{2}_{q}(\mathcal{X})$. Then by the Hilbert projection theorem (Tsiatis,
2006), there exists a unique projection $b_{\tilde{\theta}}^{*}$ of any
$b_{\tilde{\theta}}\in L^{2}_{q}(\mathcal{X})$. Define the projection as

$T_{\mathcal{F}_{\tilde{\theta}}}(b_{\tilde{\theta}})=\\{b_{\tilde{\theta}}^{*}\in\mathcal{F}_{\tilde{\theta}}:\|b_{\tilde{\theta}}^{*}-b_{\tilde{\theta}}\|_{L^{2}_{q}}=\inf_{b\in\mathcal{F}_{\tilde{\theta}}}\|b-b_{\tilde{\theta}}\|_{L^{2}_{q}}\\}.$
(8)

We now describe how $T_{{\mathcal{F}}_{\tilde{\theta}}}$ is used to define the
projection posterior. Suppose a prior distribution
$\Pi(\theta,b)=\Pi(\theta)\Pi(b)$ is elicited on $\Theta\times
L^{2}_{q}(\mathcal{X})$. Let
$\widetilde{B}=\widetilde{B}_{1}\times\widetilde{B}_{2}$ be a measurable
subset of the Borel $\sigma$-algebra of
$\Theta\times{\mathcal{F}}_{\tilde{\theta}}$ . Then given $y^{(n)}$, we define
the posterior probability $\Pi_{\text{proj}}(\widetilde{B}\mid y^{(n)})$ of
$\widetilde{B}$ under the prior $\Pi(\theta,b)$ as $\Pi(B\mid y^{(n)})$ where
$B=\widetilde{B}_{1}\times
T_{\mathcal{F}_{\tilde{\theta}}}^{-1}(\widetilde{B}_{2})$ where
$T_{\widetilde{\mathcal{F}}_{\theta}}^{-1}(\widetilde{B}_{2})=\\{b_{\tilde{\theta}}:T_{\mathcal{F}_{\tilde{\theta}}}(b_{\theta})\in\widetilde{B}_{2}\\}$,
that is,

$\Pi_{\text{proj}}(\widetilde{B}\mid y^{(n)})=\Pi(B\mid
y^{(n)})=\dfrac{\int_{B}\ell(y^{(n)}\mid\theta,b)d\Pi(\theta,b)}{\int\ell(y^{(n)}\mid\theta,b)d\Pi(\theta,b)}.$
(9)

Measurability of $T$ is guaranteed since $\mathcal{F}_{\tilde{\theta}}$ is
non-empty, closed and convex; see also Sen et al. (2018). Although (8) is
defined as a solution to an optimization problem, in this particular case, it
has an explicit form which we call functional projection.

######  Lemma 3.3.

Fix $b\in L^{2}_{q}(\mathcal{X})$ and let $\tilde{\theta}$ be defined as
above. Let the functions $g_{j,k}(x),\,j=1,\ldots,p$ satisfy the following:
$\sum_{j=1}^{p}\alpha_{j}g_{j,k}(x)=0$ for all $x\in\mathcal{X}$ and for all
$k=1,\ldots,q$ iff $\alpha_{j}=0,\,j=1,\ldots,p$. Then
$T_{\mathcal{F}_{\tilde{\theta}}}(b)=b^{*}(x)=b(x)-\sum_{j=1}^{p}\lambda_{j}^{*}g_{j}(x)$
where the vector
$\lambda=(\lambda_{1}^{*},\ldots,\lambda_{p}^{*})^{\mathrm{\scriptscriptstyle{T}}}$
satisfies $Q\lambda=\eta$ with $\eta=(\sum_{k=1}^{q}\langle
b_{k},g_{1,k}\rangle,\ldots,\sum_{k=1}^{q}\langle
b_{k},g_{p,k}\rangle)^{\mathrm{\scriptscriptstyle{T}}}$ and $Q$ is a $p\times
p$ matrix with elements $Q_{jj^{\prime}}=\sum_{k=1}^{q}\langle
g_{j,k},g_{j^{\prime},k}\rangle$.

The proof is provided in the Appendix. Having defined the projection for our
purpose, we can then devise an MCMC algorithm to sample from
$\Pi_{\text{proj}}(\cdot\mid y^{(n)})$ defined in (9).

1\. Update $b_{\tilde{\theta}}\sim\Pi(b_{\tilde{\theta}}\mid\theta,y^{(n)})$

2\. Project $b_{\tilde{\theta}}$ onto $\mathcal{F}_{\tilde{\theta}}$ following
the Proposition 3.3 to obtain $b^{*}_{\tilde{\theta}}$.

3\. Update $\theta\sim\Pi(\theta\mid b^{*}_{\tilde{\theta}},y^{(n)})$.

Algorithm 1 Projection sampler 1

As an example, consider the case when $b_{\tilde{\theta}}$ is endowed with a
zero mean Gaussian process prior with covariance kernel $C(\cdot,\cdot)$. Then
following Algorithm 1, we update
$b_{\tilde{\theta}}\sim\Pi(b\mid\tilde{\theta},y^{(n)})$ which is a
multivariate Gaussian distribution. Next, we compute the corresponding
projection $b^{*}_{\tilde{\theta}}$ using Lemma 3.3. Finally, we sample
$\theta\sim\Pi(\theta\mid b^{*}_{\tilde{\theta}},y^{(n)})$. This strategy also
works for other priors such as BART, where in the first step we update the
parameters of the BART parameters using their respective full conditionals.
The next steps are exactly the same. An alternative finite-dimensional
projection method is described in the supplementary document.

For a full Bayesian inference on (2) under the projection posterior framework
described above, one should ideally elicit a prior $\Pi(\Sigma_{F})$ on the
unknown covariance matrix $\Sigma_{F}$. However, this often adds to
computational challenges already present in a calibration problem and may not
influence the uncertainty observed in $\theta$ \- the central goal of these
problems (Bayarri et al., 2007). As an alternative, a modular Bayes approach
is considered here where we use a plug-in estimate $\hat{\Sigma}_{F}$. This
estimate is constructed by fitting a nonparametric regression model, i.e.
$y_{F,k}(x)=y_{R,k}(x)+\epsilon_{k}(x)$ to each outcome separately with
idiosyncratic variance parameters $\sigma_{F,k}^{2}$ with conjugate Inverse
Gamma priors and $y_{R,k}\sim\Pi_{k}$ for some nonparametric prior on
$y_{R,k}$. We then set $\hat{\Sigma}_{F}=\text{Cov}(E)$ where $E$ is the error
matrix with columns $E_{k}=y_{F,k}-\bar{y}_{R}$ where $\bar{y}_{R}$ is the
posterior mean of predictions on the training set. Given $\hat{\Sigma}_{F}$,
it can be assumed without loss of generality that
$\epsilon(x)\sim{\mathrm{N}}(0,\mathrm{I}_{q})$ in (2). This also makes prior
elicitation on the multivariate bias $b$ simpler; a natural choice is
$\Pi(b)=\prod_{k=1}^{q}\Pi(b_{k})$.

##  4 Explicitly unavailable computer model

Until now we have assumed that the computer model $f(x,t)$ is either
explicitly specified or can be evaluated cheaply using a code for every
$(x,t)\in\mathcal{X}\times\Theta$. However, in many calibration problems,
including ours, this is not the case; the computer code can be computationally
very expensive; see also Section 2. The standard approach in the literature
has been to approximate the computer model using a nonparametric method,
typically the Gaussian process (Santner et al., 2003). While the original
calibration framework developed by Kennedy and O’Hagan (2001) modeled the
observed data and the computer model simultaneously, a modular approach to
inference has been advocated by many authors, e.g. Bayarri et al. (2007);
Plumlee (2017) wherein the modeling of the explicitly defined computer model
is done separately from the modeling of the observed data. Here, we adopt this
modular approach with one crucial difference with Plumlee (2017) where the
author uses a probabilistic definition of the computer model using a Gaussian
process. We, on the other hand, use a deterministic definition $\hat{f}$ of
the computer model obtained using any nonparametric approximator such as the
posterior mean of a BART fit, a neural net obtained by minimizing the least
square error between code evaluations and the neural net output, or a random
forest fit; such an approach was also adopted by Xie and Xu (2021). Our
motivation for doing so is to allow for more computational and modeling
flexibility than a Gaussian process framework. Specifically, we work with the
posterior mean of a BART fit which proved to be superior among all other
choices of models for $f(x,t)$ in our numerical experiments as well as in our
motivating example. In particular, we compared the squared error loss on a
held-out test data of size $n_{t}=50$ for these methods:
$(n_{t})^{-1}\sum_{i=1}^{n_{t}}\\{f(x,t_{i})-\hat{f}(x,t_{i})\\}^{2}$ for each
location $x$ in the observed data $x\in\mathcal{X}$ and each output of the
computer model. In the best case scenario, the ratio of the squared error loss
for the posterior predictive mean improved over the second best method (Random
Forest) by a factor of almost 5.

##  5 Simulations and univariate applications

###  5.1 Univariate simulations

In this section, we compare our method to the GP orthogonal calibration by
Plumlee (2017) (OGP) and the projected calibration (PCAL) method by Xie and Xu
(2021) in several simulation experiments. Since these methods were originally
developed for one outcome, we restrict here to the case $q=1$. We focus on
each of these methods’ ability to estimate $\theta^{*}$, uncertainty
quantification, and their associated computing time. For generating the data,
we consider two situations when the definition of the computer model is
available explicitly:

  1. 1.

Model 1: $y_{R}(x)=4x+x\sin 5x$, $f(x,t)=tx$ where $x\in[0,1]$ (Plumlee, 2017,
Example 5.1),

  2. 2.

Model 2: $f(x,t)=7\\{\sin(2\pi t_{1}-\pi)\\}^{2}+2\\{(2\pi
t_{2}-\pi)^{2}\sin(2\pi x-\pi)\\}$, $y_{R}(x)=f(x,\theta^{*})$ where
$\theta^{*}=(0.2,0.3)$, $x\in[0,1]$ (Xie and Xu, 2021, Configuration 1).

In Model 1, $\theta^{*}=3.56$. We additionally consider one more situation
where we do not assume the knowledge of $f(x,t)$ in Model 2 and approximate it
by the posterior predictive mean of a BART when $N=7^{2}$ code outputs are
available on the observed values of $x$, meaning we have for each value
$\\{\theta_{1},\ldots,\theta_{N}\\}$ the computer code output is available for
$f(\theta_{k},x_{i}),\,k=1,\ldots,N,i=1,\ldots,n$.

  * •

Model 3: $f(x,t)$ is same as Model 2 but only outputs $f_{1},\ldots,f_{N}$ are
available. Everything else is same as Model 2.

In this case, we set $\tilde{\theta}$ to be the parameter value at which the
sample $L^{2}$-distance is minimized between the observed data and the
posterior predictive mean of $f$.

For all cases we simulate $n=100$ field observations from the model
$y^{F}(x_{i})=y^{R}(x_{i})+\epsilon_{i},\,\epsilon_{i}\sim{\mathrm{N}}(0,0.2^{2})$
where we sample the $x_{i}$’s uniformly over $[0,1]$. For the proposed method,
we consider two priors on the bias function - GP and BART. For the GP prior,
we use the covariance kernel
$C(x,x^{\prime})=\sigma^{2}(1+|x-x^{\prime}|/\psi)\exp\\{-|x-x^{\prime}|/\psi\\}$.
We set $\psi=1/2$ and use a plug-in estimate of $\sigma^{2}$ which is obtained
by fitting a non-parametric BART regression model to the observed field data
$(y_{F,i},x_{i})_{i=1}^{n}$. The corresponding projected versions are
abbreviated as PGP and PBART. We also consider the functional and finite
dimensional projections for PGP and PBART, with the finite-dimensional version
abbreviated as PGP(Fn) and PBART(Fn). For Model 3, we only consider the
infinite dimensional projections. In Tables 1, 2 and 3 we report the results
of our experiments from a replicated simulation study. All numbers are
averaged over 100 replications and for each model we draw 5000 posterior
samples using MCMC with 1000 burn in samples. In terms of estimation quality
of the calibration parameters, all the different methods perform similarly.
Indeed, the posterior mean and standard deviation in Table 1 and Table 2 are
almost identical. However, the key difference is in the respective runtimes of
these methods. The runtimes reported in the tables are comparable up to coding
language and efficiencies of these implementations. Nonetheless, PBART
outperforms all other methods by quite a large margin providing at least an
order of magnitude improvement over the next fastest method. We note here, the
runtimes reported for PCAL correspond to the standard implementation of the
method and not the approximated calibration method proposed in Xie and Xu
(2021); the approximated calibration method is much faster compared to its
standard version although even for that the average runtime for Model 2 is
36.33 seconds which is almost 4 times the time of PBART. In terms of
uncertainty quantification, all the methods in consideration provide roughly
the nominal coverage although it seems that PGP and PBART is biased downwards
(but only slightly). This is also apparent in the relatively sharper posterior
distributions of $\theta$ plotted in Figure 2.

| Model 1 ($\theta^{*}=3.56$)
---|---
| PGP | PGP(Fn) | PBART | PBART(Fn) | OGP | PCAL
Mean | 3.56 | 3.59 | 3.57 | 3.58 | 3.57 | 3.57
Std. Dev. | 0.008 | 0.007 | 0.02 | 0.02 | 0.003 | 0.003
Coverage | 0.93 | 0.93 | 0.92 | 0.93 | 0.95 | 0.95
Runtime | 3.78 sec | 30.72 sec | 1.86 sec | 17 min | 1.88 | 13.92 min

Table 1: Simulation results for Model 1. We report the posterior mean,
posterior standard deviation, coverage of the 95% credible intervals and the
average runtime.

| Model 2 ($\theta^{*}=(0.2,0.3)$) |
---|---|---
| PGP | PGP(Fn) | PBART | PBART(Fn) | OGP | PCAL
Mean | (0.2, 0.3) | (0.2, 0.31) | (0.19, 0.3) | (0.2, 0.31) | (0.2, 0.3) | (0.2, 0.3)
Std. Dev. | (0.0002, 0.0002) | (0.0003, 0.0004) | (0.0008, 0.0007) | (0.0008, 0.0008) | (0.001, 0.002) | (0.003, 0.0006)
Coverage | 0.91 | 0.91 | 0.92 | 0.92 | 0.96 | 0.96
Runtime | 18.18 sec | 2.65 min | 9.38 sec | 1.35 hrs | 16.58 min |

Table 2: Simulation results for Model 2. We report the posterior mean,
posterior variance, coverage of the 95% credible intervals and the average
runtime.

| Model 3 ($\theta_{0}=(0.2,0.3)$)
---|---
| PGP | PBART
Mean | (0.21, 0.33) | (0.22, 0.34)
Std. Dev. | (0.03, 0.03) | (0.04, 0.03)
Runtime | 24 min | 20 min

Table 3: Simulation results for Model 3. We report the posterior mean,
posterior variance and the average runtime.

(a) Model 1: $\pi(\theta\mid\text{data})$

(b) Model 2: $\pi(\theta_{1}\mid\text{data})$

(c) Model 2: $\pi(\theta_{2}\mid\text{data})$

Figure 2: Posterior distribution of the calibration parameter $\theta$.

###  5.2 Univariate analysis on Schwarzschild model outputs

We apply the proposed posterior projection technique on each output of the
Schwarzschild model separately. We illustrate the specifics of our application
with the first output of the Schwarzschild model, which is the mean velocity.
Similar to the Kennedy and O’Hagan (2001) setup, we have the following model
for the observed velocities $v_{F}(x)$ and the corresponding code output
$v_{S}(x)$ at location $x\in\mathcal{X}$:

$v_{F}(x_{i})=v_{R}(x_{i})+\epsilon_{i}=v_{S}(x_{i};\tilde{\theta})+b_{v,\tilde{\theta}}(x_{i})+\epsilon_{i},\quad\epsilon_{i}\sim{\mathrm{N}}(0,\sigma_{v}^{2}),$
(10)

where we let $\tilde{\theta}$ to be the parameter value where the $L^{2}$ loss
defined in previous sections is minimized and $b_{v,\tilde{\theta}}(\cdot)$ is
the bias function corresponding to mean velocity at $\tilde{\theta}$. Since
the explicit form of the model is not available, we use the posterior mean of
a BART fit as the definition of $v_{S}(x;t)$ based on the model outputs only.
We fix $\sigma_{v}$ to the posterior mean of $\sigma$ of a BART regression fit
to $(v_{F,i},x_{i})_{i=1}^{n}$. We set
$\theta_{j}\sim{\mathrm{N}}(0,\gamma^{2})$ with $\gamma=10$ as priors for
$\theta_{j}$, and $\Pi(\theta)=\prod_{j=1}^{p}\Pi(\theta_{j})$. We implement
PGP and PBART under these settings where for implementing Step 3 of Algorithm
1, we use an adaptive Metropolis-Hastings sampler Haario et al. (2001) which
is tuned to maintain an acceptance rate of approximately 0.3.

In Table 4, we summarize the posterior distribution of the calibration
parameters for all the four different outputs of the Schwarzschild model.
Specifically, we report the posterior mean, standard deviation and the 95%
credible interval for all the three calibration parameters
$(\theta_{1},\theta_{2},\theta_{3})$. As expected, there are some differences
in the posterior means for different velocity moments being fit. For example,
the posterior mean for $\theta_{1}$ differs when using $h_{4}$ compared to
when using $\tau$. The difference may be due to the physical degeneracy
between the orbital anisotropy and black hole mass, such that large $\tau$ at
the center of the 2-dimensional spatial grid can imply a large black hole mass
but so too can a smaller $\tau$ with radial anisotropy, which is reflected in
the values of $h_{4}$. Or, it could be due simply to the model inadequacy of a
univariate treatment of the problem. In addition, there seems to a broad
agreement between the posterior means for $\theta_{1}$ and $\theta_{2}$ under
the two priors. However, the key difference is in the posterior mean of
$\theta_{3}$. For the GP prior, the posterior mean for the 4 outcomes is
roughly around 1, whereas for BART, the posterior mean is around 1.4. One
possible explanation for this discrepancy is that for $\theta_{3}$, the
computer code output is only available for $\theta_{3}=1,2,3$. The scarcity of
data along this dimension may potentially lead to unstable estimates. A joint
model accounting for the covariance between the measurements can reconcile
this apparent discrepancy in the posterior distribution and the problems
induced by sparse data. Numerical experiments in the next subsection
illustrates the benefits of a joint model when outcomes are correlated.

|  | PGP | PBART
---|---|---|---
|  | $\theta_{1}$ | $\theta_{2}$ | $\theta_{3}$ | $\theta_{1}$ | $\theta_{2}$ | $\theta_{3}$
$v$ | Mean | 9.46 | 9.06 | 0.82 | 9.59 | 9.16 | 1.22
Std. Dev. | 0.14 | 0.08 | 0.09 | 0.27 | 0.15 | 0.07
Interval | (9.21, 9.80) | (8.93, 9.18) | (0.65, 1.00) | (9.20, 9.97) | (8.93, 9.50) | (1.12, 1.32)
$\tau$ | Mean | 9.98 | 9.12 | 1.02 | 9.59 | 9.03 | 1.42
Std. Dev. | 0.09 | 0.07 | 0.10 | 0.10 | 0.13 | 0.18
Interval | (9.81, 10.14) | (8.99, 9.26) | (0.87, 1.23) | (9.41, 9.76) | (8.84, 9.31) | (1.15, 1.76)
$h_{3}$ | Mean | 9.96 | 9.36 | 1.04 | 10.09 | 9.27 | 1.41
Std. Dev. | 0.07 | 0.08 | 0.12 | 0.07 | 0.08 | 0.06
Interval | (9.86, 10.11) | (9.15, 9.49) | (0.87, 1.25) | (9.97, 10.23) | (9.12, 9.41) | (1.25 , 1.41)
$h_{4}$ | Mean | 10.18 | 9.07 | 1.09 | 9.73 | 9.09 | 1.53
Std. Dev. | 0.11 | 0.06 | 0.06 | 0.06 | 0.12 | 0.09
Interval | (9.89, 10.62) | (8.97, 9.20) | (0.97, 1.19) | (9.61, 9.82) | (8.85, 9.28) | (1.35, 1.69)

Table 4: Posterior summaries for the calibration parameters $\theta_{1}$ =
logarithm of mass of the black hole, $\theta_{2}$ = stellar mass-to-light
ratio and $\theta_{3}$ = fraction of dark matter. The posterior mean, standard
deviation (Std. Dev.) and the 95% credible interval is reported here.

###  5.3 Multivariate simulations

Here we consider a simulation scenario where more than outcome is measured,
i.e. $q>1$. We illustrate this with the case $q=2$. For $x\in[0,1]$, consider
the model

$\begin{pmatrix}y_{F,1}(x)\\\
y_{F,2}(x)\end{pmatrix}=\begin{pmatrix}f_{1}(x,\theta)\\\
f_{2}(x,\theta)\end{pmatrix}+\begin{pmatrix}b_{\theta,1}(x)\\\
b_{\theta,2}(x)\end{pmatrix}+\begin{pmatrix}\epsilon_{1}\\\
\epsilon_{2}\end{pmatrix}$ (11)

where $y_{R,1}(x)=4x+x\sin 5x$, $y_{R,2}(x)=[1+\exp\\{-6(x-0.5)\\}]^{-1}$ are
the real data-generating processes. The respective computer models are
$f_{1}(x,t)=tx$, $f_{2}(x,t)=\Phi\\{t(x-0.5)\\}$ with $\Phi(\cdot)$ being the
cdf of $Z\sim{\mathrm{N}}(0,1)$. The errors are assumed to follow
${\mathrm{N}}(0,\Sigma)$. The loss
$L\\{y_{R}(\cdot),f(\cdot,t)\\}=\sum_{k=1}^{2}\int_{\mathcal{X}}\\{y_{R,k}(x)-f_{k}(x,t)\\}^{2}$
is minimized at $\theta^{*}=3.56$ approximately. However, the behavior of the
$L_{2}$-loss in these two cases is very different, especially near the
minimum. Indeed, as shown in the left panel of Figure 3, for the pair
$(y_{R,1},f_{1})$, the loss increases quite sharply around 3.56, whereas for
$(y_{R,2},f_{2})$, the increase in loss for $t>3.56$ is much slower. As such,
when data from model (11) is fitted separately in the univariate framework,
one should expect more uncertainty in $\theta$ for $(y_{R,2},f_{2})$. On the
other hand, a joint multivariate fit should reflect the fact when the two real
processes and their corresponding computer models are compared together, the
posterior distribution of $\theta$ is much more centered around 3.56.

(a) $L(\theta)=\int_{0}^{1}\left\\{y_{R,1}(x)-f_{1}(x,t)\right\\}^{2}$

(b) $L(\theta)=\int_{0}^{1}\left\\{y_{R,2}(x)-f_{2}(x,t)\right\\}^{2}$

Figure 3: Loss comparison

For data generation, we consider $n=100$ and a covariance matrix $\Sigma_{0}$
such that the diagonal elements are $0.2^{2}$ and the off-diagonal elements
are $0.012$. We consider PGP and PBART for this setup and use the functional
projection. We fit this data separately in the univariate framework and then a
joint model is fitted. For a fair comparison all hyperparameter values were
the same. We plot the posterior distribution of $\theta$ obtained using the
PBART method in Figure 4. In the first and second panel, the posterior
distribution of $\theta$ is shown when two separate univariate models are
fitted, and in the third panel we show that for a multivariate fit. As
mentioned earlier, the posterior distribution of $\theta$ when data for the
second outcome is fitted separately, is more dispersed and has relatively
large mass in the interval [3, 4]. This is not the case for the first outcome,
as almost the entire mass is within [3.5, 3.65]. However, the posterior
distribution of $\theta$ form a multivariate fit is able to borrow information
across the two outcomes so that most of the mass of the distribution is within
[3.4,3.7].

Figure 4: Posterior distribution of $\theta$ for univariate and multivariate
PBART fits for model 11. Left two panels show posteriors for two separate
univariate fits, and the third panel shows posterior from a combined
multivariate fit.

##  6 Calibrating the Schwarzschild model

Here we carry out a joint calibration of the Schwarzschild model. For this, we
first estimate the covariance matrix as $\Sigma=\text{Cov}(E)$, where $E$ is a
$n\times q$ error matrix obtained by fitting univariate BART regressions to
the field data of each outcome. Here we fit the PBART method with functional
projection. The marginal posterior distribution of $\theta$ is shown in Figure
5. The 95% symmetric posterior (marginal) credible intervals are [9.54, 9.87],
[8.78, 9.13] and [0.97, 1.38] for $\theta_{1},\theta_{2},\theta_{3}$,
respectively. Intervals for $\theta_{1},\theta_{2}$ from the joint analysis
clearly seem to lie at the intersection of the intervals from the univariate
model, but the posterior of $\theta_{3}$ exhibits slight bimodality. Overall,
by modeling the four outcomes simultaneously uncertainty in $\theta$ have
reduced which may significantly cut down computation time for future
applications with an expanded parameter space. We also show the marginal
posterior predictive mean of the outcomes. In the first panel of Figure 6, we
plot the posterior predictive mean at the observed locations for
$v_{R}(\cdot)$. For reference, the observed data is plotted on the third panel
of Figure 6. Similar plots for $\tau_{R}(\cdot),h_{3,R}(\cdot)$ and
$h_{4,R}(\cdot)$ are provided in Figures 7, 8 and 9, respectively. As
mentioned earlier, the outcomes in this context represent the first four
moments of velocity, and the model performance decreases as the moments
increase. We also observed a similar phenomenon when we looked at the relative
squared error loss for a held-out data set of size 20. For $v_{F}(x)$, We
define the relative squared error loss as
$n_{t}^{-1}\sum_{i=1}^{n_{t}}(1-\hat{v}_{F}(x_{i})/v_{F}(x_{i}))^{2}$ where
$\hat{v_{F}}(x_{i})$ represents the posterior mean of samples of
$[v_{S}(x_{i},\theta)+b_{v,\theta}(x_{i})]$. Here $n_{t}$ is the test data
size. The definition is similar for other outcomes. The values we obtained are
0.17, 0.003, 3.92 and 30.42 for $v,\tau,h_{3},h_{4}$, respectively.

Figure 5: Posterior distribution of $\theta$ from multivariate PBART fit for
the Schwarzschild model.

Figure 6: Posterior predictive mean of PBART is plotted for $v_{R}(\cdot)$ in
the first panel. In the second panel, the observed values $v_{F}(\cdot)$ are
also shown for reference.

Figure 7: Posterior predictive mean of PBART is plotted for $\tau_{R}(\cdot)$
in the first panel. In the second panel, the observed values $\tau_{F}(\cdot)$
are also shown for reference.

Figure 8: Posterior predictive mean of PBART is plotted for $h_{3,R}(\cdot)$
in the first panel. In the second panel, the observed values $h_{3,F}(\cdot)$
are also shown for reference.

Figure 9: Posterior predictive mean of PBART is plotted for $h_{4,R}(\cdot)$
in the first panel. In the second panel, the observed values $h_{4,R}(\cdot)$
are also shown for reference.

##  7 Discussion

The Schwarzschild model, although computationally intensive, is an important
tool in understanding the dynamics of a black hole and its host galaxy.
Motivated by this application we developed a multivariate calibration method
that ensures parameter identifiablity under the squared error loss. The key
benefits of the proposed projection posterior approach is its flexibility in
accommodating user-specified prior distributions on the multivariate bias
function, and the fact that such a projection is available analytically.
Benefits of a multivariate analysis is demonstrated through numerical
experiments and it seems to impact the analysis of Schwarzschild model
positively as well. Conceptually, the projection approach can be extended to
alternative loss functions that lead to different sufficiency conditions for
identifiability and high-dimensional outcomes where low-rank models for the
error covariance might be warranted.

## References

  * Bayarri et al. (2007) Bayarri, M., D. Walsh, J. Berger, J. Cafeo, G. Garcia-Donato, F. Liu, J. Palomo, R. Parthasarathy, R. Paulo, and J. Sacks (2007).  Computer model validation with functional output.  Annals of Statistics 35(5), 1874–1906.
  * Bayarri et al. (2007) Bayarri, M. J., J. O. Berger, R. Paulo, J. Sacks, J. A. Cafeo, J. Cavendish, C.-H. Lin, and J. Tu (2007).  A framework for validation of computer models.  Technometrics 49(2), 138–154.
  * Chakraborty et al. (2013) Chakraborty, A., B. K. Mallick, R. G. Mcclarren, C. C. Kuranz, D. Bingham, M. J. Grosskopf, E. M. Rutter, H. F. Stripling, and R. P. Drake (2013).  Spline-based emulators for radiative shock experiments with measurement error.  Journal of the American Statistical Association 108(502), 411–428.
  * Chipman et al. (2010) Chipman, H. A., E. I. George, and R. E. McCulloch (2010).  BART: Bayesian additive regression trees.  The Annals of Applied Statistics 4(1), 266–298.
  * Denison et al. (1998) Denison, D. G., B. K. Mallick, and A. F. Smith (1998).  Bayesian mars.  Statistics and Computing 8(4), 337–346.
  * Forest et al. (2008) Forest, C. E., B. Sansó, and D. Zantedeschi (2008).  Inferring climate system properties using a computer model.  Bayesian Analysis 3(1), 1–37.
  * Friedman (1991) Friedman, J. H. (1991).  Multivariate adaptive regression splines.  The annals of statistics 19(1), 1–67.
  * Gattiker et al. (2006) Gattiker, J., D. Higdon, S. Keller-McNulty, M. McKay, L. Moore, and B. Williams (2006).  Combining experimental data and computer simulations, with an application to flyer plate experiments.  Bayesian Analysis 1(4), 765–792.
  * Gneiting and Raftery (2007) Gneiting, T. and A. E. Raftery (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American statistical Association 102(477), 359–378.
  * Haario et al. (2001) Haario, H., E. Saksman, J. Tamminen, et al. (2001).  An adaptive Metropolis algorithm.  Bernoulli 7(2), 223–242.
  * Higdon et al. (2008) Higdon, D., J. Gattiker, B. Williams, and M. Rightley (2008).  Computer model calibration using high-dimensional output.  Journal of the American Statistical Association 103(482), 570–583.
  * Kennedy and O’Hagan (2001) Kennedy, M. C. and A. O’Hagan (2001).  Bayesian calibration of computer models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63(3), 425–464.
  * Kormendy and Ho (2013) Kormendy, J. and L. C. Ho (2013, August).  Coevolution (Or Not) of Supermassive Black Holes and Host Galaxies.  Annual Review of Astronomy and Astrophysics 51(1), 511–653.
  * Linero and Yang (2018) Linero, A. R. and Y. Yang (2018).  Bayesian regression tree ensembles that adapt to smoothness and sparsity.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 80(5), 1087–1110.
  * Marmin and Filippone (2022) Marmin, S. and M. Filippone (2022).  Deep gaussian processes for calibration of computer models.  Bayesian Analysis 1(1), 1–30.
  * McConnell and Ma (2013) McConnell, N. J. and C.-P. Ma (2013, February).  Revisiting the Scaling Relations of Black Hole Masses and Host Galaxy Properties.  The Astrophysical Journal 764(2), 184.
  * Mehrgan et al. (2019) Mehrgan, K., J. Thomas, R. Saglia, X. Mazzalay, P. Erwin, R. Bender, M. Kluge, and M. Fabricius (2019, December).  A 40 Billion Solar-mass Black Hole in the Extreme Core of Holm 15A, the Central Galaxy of Abell 85.  The Astrophysical Journal 887(2), 195.
  * Neal (2012) Neal, R. M. (2012).  Bayesian learning for neural networks, Volume 118.  Springer Science & Business Media.
  * Plumlee (2017) Plumlee, M. (2017).  Bayesian calibration of inexact computer models.  Journal of the American Statistical Association 112(519), 1274–1285.
  * Plumlee and Joseph (2018) Plumlee, M. and V. R. Joseph (2018).  Orthogonal gaussian process models.  Statistica Sinica, 601–619.
  * Pollard and Radchenko (2006) Pollard, D. and P. Radchenko (2006).  Nonlinear least-squares estimation.  Journal of Multivariate Analysis 97(2), 548–562.
  * Pratola and Higdon (2016) Pratola, M. T. and D. M. Higdon (2016).  Bayesian additive regression tree calibration of complex high-dimensional computer models.  Technometrics 58(2), 166–179.
  * Quenneville et al. (2021) Quenneville, M. E., C. M. Liepold, and C.-P. Ma (2021).  Dynamical modeling of galaxies and supermassive black holes: Axisymmetry in triaxial schwarzschild orbit superposition models.  The Astrophysical Journal Supplement Series 254(2), 25\.
  * Rasmussen and Williams (2005) Rasmussen, C. E. and C. K. I. Williams (2005).  Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning).  The MIT Press.
  * Ročková and van der Pas (2020) Ročková, V. and S. van der Pas (2020).  Posterior concentration for bayesian regression trees and forests.  The Annals of Statistics 48(4), 2108–2131.
  * Sacks et al. (1989) Sacks, J., W. J. Welch, T. J. Mitchell, and H. P. Wynn (1989).  Design and analysis of computer experiments.  Statistical science 4(4), 409–423.
  * Saglia et al. (2016) Saglia, R. P., M. Opitsch, P. Erwin, J. Thomas, A. Beifiori, M. Fabricius, X. Mazzalay, N. Nowak, S. P. Rusli, and R. Bender (2016, February).  The SINFONI Black Hole Survey: The Black Hole Fundamental Plane Revisited and the Paths of (Co)evolution of Supermassive Black Holes and Bulges.  The Astrophysical Journal 818(1), 47.
  * Salter et al. (2018) Salter, J. M., D. B. Williamson, J. Scinocca, and V. Kharin (2018).  Uncertainty quantification for spatio-temporal computer models with calibration-optimal bases.  arXiv preprint arXiv:1801.08184 12.
  * Santner et al. (2003) Santner, T. J., B. J. Williams, and W. I. Notz (2003).  The Design and Analysis of Computer Experiments.  Springer series in statistics. Springer.
  * Schwarzschild (1979) Schwarzschild, M. (1979).  A numerical model for a triaxial stellar system in dynamical equilibrium.  The Astrophysical Journal 232, 236–247.
  * Sen et al. (2018) Sen, D., S. Patra, and D. Dunson (2018).  Constrained inference through posterior projections.  arXiv preprint arXiv:1812.05741.
  * Tsiatis (2006) Tsiatis, A. A. (2006).  Semiparametric theory and missing data, Volume 4.  Springer.
  * Tuo and Wu (2015) Tuo, R. and C. J. Wu (2015).  Efficient calibration for imperfect computer models.  Annals of Statistics 43(6), 2331–2352.
  * Tuo and Wu (2016) Tuo, R. and C. J. Wu (2016).  A theoretical framework for calibration in computer models: parametrization, estimation and convergence properties.  SIAM/ASA Journal on Uncertainty Quantification 4(1), 767–795.
  * van den Bosch et al. (2008) van den Bosch, R., G. Van De Ven, E. Verolme, M. Cappellari, and P. De Zeeuw (2008).  Triaxial orbit based galaxy models with an application to the (apparent) decoupled core galaxy ngc 4365.  Monthly Notices of the Royal Astronomical Society 385(2), 647–666.
  * Xie and Xu (2021) Xie, F. and Y. Xu (2021).  Bayesian projected calibration of computer models.  Journal of the American Statistical Association 116(536), 1965–1982.
  * Yang et al. (2017) Yang, Y., A. Bhattacharya, and D. Pati (2017).  Frequentist convergence and sup-norm convergence rate in gaussian process regression.  arxiv preprint arXiv:1708.04753.

##  Appendix A Proofs

###  A.1 Proof of Proposition 3.2

###### Proof.

We first characterize $\theta$ in terms of the expected loss. Recall that for
$Q=\mathrm{I}_{q}$,
$\theta=\operatorname*{arg\,min}_{t\in\Theta}\int_{x}\sum_{k=1}^{q}\\{y_{R,k}(x)-f_{k}(x,t)\\}^{2}dx$.
For any fixed $t\in\Theta$, let $e(x)=\\{y^{F}(x)-f(x,t)\\}$, then

$\displaystyle
E_{P}\\{e(x)^{\mathrm{\scriptscriptstyle{T}}}e(x)\\}=(\text{Vol}(\mathcal{X}))^{-1}\int_{x}\left[\sum_{k=1}^{q}E_{P_{y\mid
x}}\\{e^{2}_{k}(x)\\}\right]dx$
$\displaystyle=(\text{Vol}(\mathcal{X}))^{-1}\int_{x}\sum_{k=1}^{q}E_{P_{y\mid
x}}\\{y_{F,k}-y_{R,k}(x)\\}^{2}dx+(\text{Vol}(\mathcal{X}))^{-1}\int_{x}\sum_{k=1}^{q}\\{y_{R,k}(x)-f_{k}(x,t)\\}^{2}dx$
$\displaystyle=\text{tr}(\Sigma_{F})+(\text{Vol}(\mathcal{X}))^{-1}\int_{x}e(x)^{\mathrm{\scriptscriptstyle{T}}}e(x)dx,$

where $\text{tr}(A)$ denotes the trace of a matrix $A$. Hence, $\theta$ can be
seen as

$\operatorname*{arg\,min}_{t\in\Theta}\int_{x}\sum_{k=1}^{q}\\{y_{R,k}(x)-f_{k}(x,t)\\}^{2}dx=\operatorname*{arg\,min}_{t\in\Theta}E_{P}[e(x)^{\mathrm{\scriptscriptstyle{T}}}e(x)].$

Now define
$Q_{n,k}(t)=n^{-1}\sum_{i=1}^{n}\\{y_{F,k}(x_{i})-f_{k}(x_{i},t)\\}^{2}$ and
$Q_{n}(t)=q^{-1}\sum_{k=1}^{q}Q_{n,k}(t)$ which is simply
$E_{P_{n}}[e(x)^{\mathrm{\scriptscriptstyle{T}}}e(x)]$ where $P_{n}$ is the
corresponding joint empirical measure. Invoking Theorem 3 of Pollard and
Radchenko (2006) and the compactness of $\Theta$, we get that
$\sup_{f}|Q_{n}(t)-E_{P}(Q_{n}(t))|\overset{P}{\to}0$ which implies that
$\theta_{n}^{*}\overset{P}{\to}\theta$. ∎

###  A.2 Proof of Lemma 3.3

###### Proof.

The projection operator is the solution to the optimization problem

$\displaystyle\min\|b_{\theta}-b_{\theta}^{*}\|_{L^{2}_{q}}\quad\text{subject
to }\sum_{k=1}^{q}\langle g_{j,k},b_{\theta,k}^{*}\rangle=0,\,j=1,\ldots,p.$

Hence, the Lagrangian dual of the optimization problem is

$\displaystyle\min\sum_{k=1}^{q}\|b_{\theta,k}-b_{\theta,k}^{*}\|_{L}^{2}+\sum_{j=1}^{p}\lambda_{j}\left[\sum_{k=1}^{q}\langle
g_{j,k},b_{\theta,k}^{*}\rangle\right].$

Solving for $b_{\theta,k}^{*}$ gives
$b_{\theta,k}^{*}=b_{\theta,k}+\sum_{j=1}^{p}\lambda_{j}g_{j,k}$. Thus,
$\|b_{\theta,k}-b_{\theta,k}^{*}\|_{L^{2}_{q}}=\lambda^{\mathrm{\scriptscriptstyle{T}}}(\sum_{k=1}^{q}Q_{k})\lambda$
where
$\lambda=(\lambda_{1},\ldots,\lambda_{p})^{\mathrm{\scriptscriptstyle{T}}}$
and $Q_{k}^{p\times p}$ is a matrix with $(j,j^{\prime})$-th element $\langle
g_{j,k},g_{j^{\prime},k}\rangle$. Also, $\sum_{k=1}^{q}\langle
g_{j,k},b_{\theta,k}^{*}\rangle=\sum_{k=1}^{q}\langle
g_{j,k},b_{\theta,k}\rangle+\sum_{j^{\prime}=1}^{p}\lambda_{j^{\prime}}Q_{k,jj^{\prime}}$.
This implies that $Q\lambda=-\eta$ where $\eta\in\mathbb{R}^{p\times 1}$ with
$\eta_{j}=\sum_{k=1}^{q}\langle g_{j,k},b_{\theta,k}\rangle$ which proves the
result. ∎

Supplementary materials for “Orthogonal calibration via posterior projections
with applications to the Schwarzschild model”

##  Appendix S.1 Finite-dimensional projection

When the central focus of analysis is assessing uncertainty in $\theta$, and
$\Pi(b_{k})$ is a Gaussian process, then the projection posterior approach can
be implemented within a finite dimensional setting by ensuring
$\sum_{k=1}^{q}\mathbf{b}_{k}^{\mathrm{\scriptscriptstyle{T}}}\mathbf{g}_{j,k}=0$
where
$\mathbf{b}_{k}=(b_{k}(x_{1}),\ldots,b_{k}(x_{n}))^{\mathrm{\scriptscriptstyle{T}}}$
and $\mathbf{g}_{j,k}=(g_{j,k}(x_{1}),\ldots,g_{j,k}(x_{n}))$. This is a
consequence of posteriors of $\mathbf{b}_{k}$ being Gaussian under a Gaussian
process prior and projections of multivariate Gaussian random vectors to
linear subspaces. For example, suppose $X^{d\times
1}\sim{\mathrm{N}}(\mu,\Sigma)$ and consider the set
$S=\\{x:A^{\mathrm{\scriptscriptstyle{T}}}x=0\\}\subset\mathbb{R}^{n}$ where
$A^{d\times p}$ is matrix of rank $p$. The conditional distribution of $X\mid
X\in S$ is again a multivariate Gaussian with mean $\mu-\Sigma
A^{\mathrm{\scriptscriptstyle{T}}}(A\Sigma
A^{\mathrm{\scriptscriptstyle{T}}})^{-1}A\mu$ and covariance matrix
$\Sigma-\Sigma A^{\mathrm{\scriptscriptstyle{T}}}(A\Sigma
A^{\mathrm{\scriptscriptstyle{T}}})^{-1}A\Sigma$. Next, define the following
projection for any $x\in\mathbb{R}^{d}$

$P_{S}(x)=\Sigma^{1/2}\operatorname*{arg\,min}_{y\in S}\|y-\Sigma^{-1/2}x\|.$
(S.1)

Standard linear algebra yields that $P_{S}(x)$ has the form
$\Sigma^{1/2}(I-P_{A})\Sigma^{-1/2}y$ where
$P_{A}=A(A^{\mathrm{\scriptscriptstyle{T}}}A)^{-1}A^{\mathrm{\scriptscriptstyle{T}}}$
is the projection matrix of $A$. This leads to the following interpretation of
the multivariate Gaussian distribution conditional on linear equality
constraints.

######  Proposition S.1.1.

Suppose $X\sim{\mathrm{N}}(\mu,\Sigma)$ and consider the set
$S=\\{x:A^{\mathrm{\scriptscriptstyle{T}}}x=0\\}\subset\mathbb{R}^{d}$ where
$A^{d\times p}$ is matrix of rank $p$. Then the random variable

$P_{S}(X)\sim{\mathrm{N}}(\mu-\Sigma
A^{\mathrm{\scriptscriptstyle{T}}}(A\Sigma
A^{\mathrm{\scriptscriptstyle{T}}})^{-1}A\mu,\Sigma-\Sigma
A^{\mathrm{\scriptscriptstyle{T}}}(A\Sigma
A^{\mathrm{\scriptscriptstyle{T}}})^{-1}A\Sigma).$

###### Proof.

We have

$\displaystyle\text{Cov}(X\mid X\in S)$ $\displaystyle=\Sigma-\Sigma
A^{\mathrm{\scriptscriptstyle{T}}}(A\Sigma
A^{\mathrm{\scriptscriptstyle{T}}})^{-1}A\Sigma$
$\displaystyle=\Sigma^{1/2}(I-\Sigma^{1/2}A^{\mathrm{\scriptscriptstyle{T}}}(A\Sigma
A^{\mathrm{\scriptscriptstyle{T}}})A\Sigma^{1/2})\Sigma^{1/2}$
$\displaystyle=\Sigma^{1/2}(I-P_{B})\Sigma^{1/2}=\Sigma^{1/2}(I-P_{B})^{2}\Sigma^{1/2}=\text{Cov}(P_{S}(X)),$

where $B=A\Sigma^{1/2}$ and $P_{B}=P_{A}$ since $A$ is full rank and $\Sigma$
is non-singular by assumption. A similar analysis shows that the mean also
agrees. ∎

The proof is provided in the appendix. We illustrate this within the
orthogonal calibration context. For the sake of simplicity, consider $q=1$. We
set $A^{n\times p}=(a_{1},\ldots,a_{p})$ where
$a_{j}=(g_{j}(x_{1},\tilde{\theta}),\ldots,g_{j}(x_{n},\tilde{\theta}))^{\mathrm{\scriptscriptstyle{T}}}$.
Let $b=(b(x_{1}),\ldots,b(x_{n}))^{\mathrm{\scriptscriptstyle{T}}}$ denote the
vector of bias function evaluated at $x_{1},\ldots,x_{n}$. Then the finite-
dimensional analogue of the orthogonality condition is
$A^{\mathrm{\scriptscriptstyle{T}}}b=0$. When $b_{\tilde{\theta}}(\cdot)$ is
endowed with a zero mean Gaussian process prior with covariance kernel
$C(\cdot,\cdot)$, then apriori the vector
$b_{\tilde{\theta}}\sim{\mathrm{N}}(0,K)$ with $K_{ij}=C(x_{i},x_{j}),\,1\leq
i,j\leq n$. Also,
$b_{\tilde{\theta}}\mid\theta,y^{(n)}\sim{\mathrm{N}}((K+\sigma^{2}_{F})^{-1}z^{(n)},(K+\sigma^{2}_{F})^{-1})$
where
$z^{(n)}=(y(x_{1})-f(x_{1},\theta),\ldots,y(x_{2})-f(x_{n},\theta))^{\mathrm{\scriptscriptstyle{T}}}$.
Using Proposition S.1.1, we can then obtain $b^{*}$. Hence, for Gaussian
process priors, only the second step in Algorithm 1 from the main draft will
change under the finite-dimensional projection setup.

We now demonstrate how the finite-dimensional projection can be executed when
$b_{\tilde{\theta}}$ is specified a non-Gaussian process prior. Here the
posterior of $\mathbf{b}_{k}$ is also non-Gaussian. Suppose the prior is
parameterized by $\eta$. For instance, with a BART prior, $\eta$ contains the
trees and leaf node parameters. A prior on $b_{\tilde{\theta}}$ is specified
by a prior distribution on $\eta$ which then produces a full conditional
distribution of $\eta\mid\theta,y^{(n)}$. A posterior sample of
$b_{\tilde{\theta}}\in\mathbb{R}^{n}$ can be drawn by sampling from
$\eta\mid\theta,y^{(n)}$. For many priors, the joint distribution can be well
approximated by a multivariate Gaussian with some mean and covariance
$(\beta,\Phi)$ (Yang et al., 2017). To find the approximating mean and
covariance, we simply draw $M$ independent draws of $b_{\tilde{\theta}}$ by
sampling $\\{\eta_{1},\ldots,\eta_{M}\\}$ independently from
$\eta\mid\theta,y^{(n)}$ and set
$\beta=M^{-1}\sum_{m=1}^{M}b_{\tilde{\theta},m}$ and
$\phi=(M-1)^{-1}(b_{\tilde{\theta},m}-\beta)(b_{\tilde{\theta},m}-\beta)^{\mathrm{\scriptscriptstyle{T}}}$
where $b_{\tilde{\theta},m}$ is the draw of $b_{\tilde{\theta}}$ corresponding
to $\eta_{m},\,m=1,\ldots,M$. We note here, that $M$ can be made arbitrarily
large compared to the dimension $n$ of $b_{\tilde{\theta}}$ so that the sample
mean and covariance provide reasonable approximations to their population
counterparts. This technique is summarized in the following algorithm.

1\. Sample $M$ independent copies of parameters $\eta\mid\theta,y^{(n)}$ and
consider the corresponding draws of $b_{\tilde{\theta}}$ as
$\\{b_{\tilde{\theta},1},\ldots,b_{\tilde{\theta,M}}\\}$.

2\. Compute the quantities $\beta=M^{-1}\sum_{m=1}^{M}b_{\tilde{\theta},m}$
and
$\Phi=(M-1)^{-1}(b_{\tilde{\theta},m}-\beta)(b_{\tilde{\theta},m}-\beta)^{\mathrm{\scriptscriptstyle{T}}}$.

2\. Set the projection as
$b_{\tilde{\theta}}^{*}=\Phi^{1/2}(I-P_{G})\Phi^{-1/2}b_{\tilde{\theta}}$.

3\. Update $\theta\sim\Pi(\theta\mid b^{*}_{\tilde{\theta}},y^{(n)})$.

Algorithm 2 Projection sampler 2: non-GP priors

Clearly, for non-Gaussian priors for $b_{k}$, the finite-dimensional
projection approach is computationally much more expensive. However, if the
prior is a GP, this approach is relatively less expensive to implement and
yields similar results to the functional projection approach described in the
main draft.