# Confidence Intervals on Multivariate Normal Quantiles for Environmental
Specification Development in Multi-axis Shock and Vibration Testing

Adam Wattsa, Thomas Thompsona, and Dustin Harveya CONTACT Adam Watts. Email:
<EMAIL_ADDRESS>aE-14 Test Engineering, Los Alamos National Laboratory, P.O.
Box 1663, Los Alamos, 87545, New Mexico, USA.

###### Abstract

This article describes two Monte Carlo methods for calculating confidence
intervals on cumulative density function (CDF) based multivariate normal
quantiles that allows for controlling the tail regions of a multivariate
distribution where one is most concerned about extreme responses. The CDF
based multivariate normal quantiles associated with bivariate distributions
are represented as contours and for trivariate distributions represented as
iso-surfaces. We first provide a novel methodology for an inverse problem,
characterizing the uncertainty on the $\tau^{\mathrm{th}}$ multivariate
quantile probability, when using concurrent univariate quantile probabilities.
The uncertainty on the $\tau^{\mathrm{th}}$ multivariate quantile probability
demonstrates inadequacy in univariate methods which neglect correlation
between multiple variates. Limitations of traditional multivariate normal
tolerance regions and simultaneous univariate tolerance methods are discussed
thereby necessitating the need for confidence intervals on CDF based
multivariate normal quantiles. Two Monte Carlo methods are discussed; the
first calculates the CDF over a tessellated domain followed by taking a
bootstrap confidence interval over the tessellated CDF. The CDF based
multivariate quantiles are then estimated from the CDF confidence intervals.
For the second method, only the point associated with highest probability
density along the CDF based quantile is calculated, which greatly improves the
computational speed compared to the first method. Monte Carlo simulation
studies are used to assess the performance of the various methods. Finally,
real data analysis is performed to illustrate a workflow for CDF based
multivariate normal quantiles in the domain of mechanical shock and vibration
to specify a minimum conservative test level for environmental specification.

###### keywords:

Quantile; Tolerance; Multivariate; Coverage; Multi-axis; Specification

##  1 Introduction

Statistical quality control has many facets, one of which are defining the
limits or bounds that a manufactured product must adhere to in order for that
product to meet quality criteria. In the domain of mechanical shock and
vibration, deriving the upper limit on the vibration or shock environmental
levels that could potentially result in damage to a system or component are of
concern. In single-axis testing, the test item is exposed to a single-axis
excitation and then is rotated and re-mounted to test the subsequent
orthogonal axes. Traditional shock and vibration testing only considers the
univariate case for a single-axis response variable. Multi-axis testing is
becoming more prominent as it provides excitation from more than one axes and
doesn’t require for the test body to be rotated and re-mounted. Tolerance
intervals have been utilized to specify a conservative upper bound as the
minimum testing levels for single-axis environmental tests as they can be
estimated from using sample sizes as low as $n=3$, but do not take into
account correlation between variates as they are a univariate statistic
method. The methodology presented here offers a new statistical method for
deriving these limits when working with multiple variates which are either
response locations and/or response directions while including the correlations
between all the variates.

Traditional univariate tolerance intervals were first introduced in the early
1940s by the statistician Abraham Wald(Wald and Wolfowitz 1946; Wald 1943,
1942). A normal tolerance interval is a statistical interval that contains at
least a proportion, $\beta$, of the distribution or population for a given
confidence, $\gamma$, where $\gamma=1-\alpha$ and $\alpha$ is the significance
level used to compute the confidence level. The term normal denotes that the
data can be represented as originating from a population that is normally
distributed. We will drop the normal prefix from tolerance interval and
tolerance bound for brevity. One-sided tolerance bounds are equivalent to a
one-sided confidence interval for a normal distribution quantile (Meeker et
al. 2023). The easiest example of this equivalency occurs for an upper
statistical limit. In this case, the one-sided upper tolerance bound for the
proportion of the population, $\beta$, and the one-sided confidence bound on
the $\tau^{\mathrm{th}}$ quantile probability, both imply encompassing the
population starting at $-\infty$ and ending at the same value. Figure 1
demonstrates the difference between a one-sided statistical bound and two-
sided statistical interval covering 90% of the population. One should note how
the one-sided upper bound includes the lower tail of the distribution; a one-
sided lower bound includes the upper tail of the distribution. When the
distribution’s parameters are known, a one-sided tolerance bound, or
equivalently, a one-sided confidence bound on the distribution’s quantile
associated with the $\tau^{\mathrm{th}}$ quantile probability, is just the
inverse cumulative distribution function (ICDF). The ICDF is also known as the
quantile function, denoted as $Q(\bullet)$. A tolerance interval or bounds
accounts for the uncertainty in the parameters of the population’s
distribution as it incorporates the sample size, $n$, directly into its
calculation. As the sample size approaches infinity, $n\rightarrow\infty$, a
one-sided tolerance bound converges to the ICDF, regardless of the desired
confidence due to having zero uncertainty on the parameters of the population.

Figure 1: Comparing a two-sided and one-sided statistical interval on the
bounds covering 90% of the population using the standard Normal distribution.
As the sample size approaches infinity, $n\rightarrow\infty$, the one-sided
tolerance bounds for a coverage value of $\beta=0.9$ and the one-sided
confidence interval for the quantile probability of $\tau=0.90$ both converge
to a value of $x=1.28$, see (2) and (3).

Recent simulation work from Chiang and Hsiao (2021) have confirmed the
theoretical properties of tolerance intervals. Univariate tolerance intervals
can be found among the industrial and medical domains (Kenett and Zacks 2021;
Francq et al. 2020). There still exists much confusion over the variety of
statistical intervals as many engineers are familiar with confidence intervals
only in the context of the distribution mean. However, confidence intervals
can be applied to other parameter or characteristic of a distribution that
represents a population or process. For example a confidence interval for the
distribution median, standard deviation, or quantile. In other contexts, an
enclosure interval such as a two-side tolerance interval that contain, or
envelopes, at least a specified proportion of a distribution is required.
Perhaps some of the mistakes around selecting the correct statistical interval
is because, “… statistics textbooks typically focus on the common confidence
intervals, occasionally make reference to tolerance intervals, and consider
prediction intervals only in the context of regression analysis. This is
unfortunate because in applications, tolerance intervals, prediction
intervals, and confidence intervals on distribution quantiles and on
exceedance probabilities are needed almost as frequently as the better-known
confidence intervals.” (Meeker et al. 2023).

Two-sided statistical intervals have both a finite lower endpoint and finite
upper endpoint. Each endpoint has an equal degree of uncertainty associated
with the parameter, or quantity of interest, located at each end of the
interval’s endpoints (Meeker et al. 2023). One-sided statistical bounds on the
other hand have a degree of uncertainty only associated with either the lower
endpoint for a one-sided lower bound or the upper endpoint for a one-sided
upper bound. One-sided statistical bounds do not have finite endpoints with a
one-sided lower bound having an upper endpoint of $+\infty$ and a one-sided
upper bound having a lower endpoint of $-\infty$. One can approximately make a
one-sided lower and one-sided upper confidence interval from a two-sided
confidence interval. For example, a 90% two-sided confidence interval is
approximately equivalent to a 95%, one-sided, upper and lower confidence
interval respectively (Meeker et al. 2023).

For univariate normal distributions, the quantile, is the value associated to
the quantile function on the $\tau^{\mathrm{th}}$ quantile probability. In
univariate statistics, the quantiles are finite scalars,
$Q(\tau)\rightarrow(-\infty,\infty):0<\tau<1$. In the multivariate case,
quantiles, are the set of vectors from the quantile function for the
$\tau^{\mathrm{th}}$ quantile probability, $Q(\tau)$. The quantiles in the
multivariate case are now an infinite set of vectors, $\boldsymbol{x}$:

$\displaystyle
Q(\tau)=\left\\{\boldsymbol{x}=(x_{1},\ldots,x_{q})\in\mathbb{R}^{q}\ :\
\tau=F(\boldsymbol{x})\right\\},$ (1)

where $F(\boldsymbol{x})$ is the CDF function for a $q$-variate normal
distribution. The term univariate quantile refer to the value associated to
the quantile function on the $\tau^{\mathrm{th}}$ quantile probability, i.e.,
$x=Q(\tau)$. The term CDF based multivariate quantile, or just CDF based
quantile, refer to the infinite set of vectors associated to the quantile
function for the $\tau^{\mathrm{th}}$ quantile probability, (1).

One-sided (upper/lower) and two-sided statistical bounds are unambiguous in
univariate statistical analysis but not in the multivariate case. We shall see
in section 3.1, traditional normal tolerance regions bound the population data
from the distribution’s center (mean) outwards (Lucagbo and Mathew 2023; Young
and Mathew 2020). These multivariate regions can be thought analogously as a
multivariate extension of the “two-sided” tolerance interval utilized in
univariate statistics. Nonetheless, there does not appear to be a multivariate
extension of the “one-sided” tolerance bound or equivalently a methodology for
constructing confidence intervals on CDF based multivariate quantiles which
would allow for controlling the tail regions of a multivariate normal
distributions. It is at these tail regions where one is most concerned about
controlling extreme responses for a given confidence value.

Tolerance regions associated with multivariate normal distributions are
thoroughly described by Chew (1966) with the most recent and canonical being
from Krishnamoorthy and Mathew (2009). Bornkamp (2018) currently provides one
of the few methodologies for calculating CDF based multivariate quantiles. The
caveat with the method is that it only produces a single, equicoordinate
vector, within the infinite set of vectors that comprise a multivariate
quantile. When used appropriately by standardizing a multivariate
distribution, the equicoordinate vector can be used very efficiently. Chen and
Welsh (2002) provide methodologies for multivariate quantiles mostly
applicable for the bivariate case but provides little in how to characterize
uncertainty on the generated quantiles. Klein and Kneib (2020) provide a
detailed methodology for CDF based bivariate quantiles by transforming the
bivariate data to the unit square but also don’t address uncertainty.

The purpose and scope of this paper is extended over several sections. Section
2 demonstrates the inadequacy of using concurrent univariate tolerance
intervals or concurrent univariate quantiles in a multivariate space. We
present a novel methodology for characterizing the uncertainty on the
$\tau^{\mathrm{th}}$ multivariate quantile probability when utilizing
concurrent univariate tolerance intervals or univariate quantiles. Section 3
presents classical tolerance regions and introduces the concept of CDF based
multivariate quantiles as either contours, iso-surfaces, or hyper iso
surfaces. Section 4 discusses the importance of standardization in
multivariate statistics. We demonstrate how to construct confidence intervals
on CDF based multivariate normal quantiles with an algorithm for the
trivariate case. The concept of a critical point, the point residing on a
multivariate quantile that has the highest relative likelihood (maximum PDF
value) of occurring, is introduced. Another algorithm is provided for
constructing confidence intervals for just the critical points on CDF based
multivariate quantiles. Confidence intervals on a critical point greatly
speeds up the computation time compared to the previous method which estimates
portions of the CDF based multivariate quantile. Monte Carlo simulation
studies and a case study in the field of mechanical shock is carried out in
section 5. Finally, we describe further research directions summarize this
paper in section 6.

Since multivariate normal distributions are the ubiquitous multivariate
distribution, we believe that developing confidence intervals on CDF based
multivariate normal quantiles and confidence intervals for critical points on
these multivariate quantiles greatly benefits a large audience, particularity
those working in field of quality control.

##  2 Inadequacy of Univariate Methods in a Multivariate Space

Univariate, one-sided, upper tolerance bounds are commonly used in quality
control. The tolerance bound is denoted below,

$\displaystyle\tilde{T}_{\beta}(1-\alpha;\beta,n)\ ,$ (2)

where $1-\alpha$ is the confidence value, $\beta$ is the proportion of the
population, and $n$ denotes the sample size. When the tilde is underneath the
variable such as, $\underaccent{\tilde}{T}$, it denotes a one-sided lower
bound. When the tilde is above the variable it denotes a one-sided, upper,
bound, such as $\tilde{T}$. Equation (2) is equivalent to a one-sided, upper,
confidence bound for the $\tau^{\mathrm{th}}$ quantile probability which is
denoted below

$\displaystyle\tilde{x}_{\tau}(1-\alpha;\tau,n)\ .$ (3)

For the example provided below, each variate corresponds to one of the three
spatial directions: X, Y, and Z. One might be inclined to use multiple, one-
sided, confidence bounds for the $\tau^{\mathrm{th}}$ quantile,
$Q(\tau_{i}):i\in\left\\{X,Y,Z\right\\}$, when working with multiple variates.
Using multiple, one-sided bounds, for each spatial direction, leads to
uncertainty as to what is the exact multivariate quantile probability,
$\tau_{{\scriptscriptstyle XYZ}}$, over the entire space because the
correlations between the variates are neglected. Looking at this predicament
through the lenses of a multiple hypothesis problem, one can bound the
multivariate quantile probability for all three variates.

For example, if $\tau_{i}=0.90:i\in\left\\{X,Y,Z\right\\}$, then the lower and
upper bounds on the multivariate quantile probability are,
$0.70\leq\tau_{{\scriptscriptstyle XYZ}}\leq 0.90$. The uncorrelated case for
this example would have a multivariate quantile probability of
$\tau_{{\scriptscriptstyle XYZ}}=(0.90)^{3}=0.729$, which is between the
bounds shown above. See appendix A for the equations to calculate these
probabilities.

The top plot of Figure 2 depicts how the multivariate quantile probability
decreases as the number of variates increases from 1 to 20. One of the most
common ways to alleviate the drop in the multivariate quantile probability is
to adjust the individual quantile probability values based on either the
Bonferroni inequality or assuming that the variates are independent, thus
uncorrelated. The bottom plot of Figure 2 depicts how the individual quantile
probability values must increase as the number of variates increases from 1 to
20 such that the multivariate quantile probability maintains a constant value,
$\tau_{\mathrm{J}}=0.90$. One can see how the individual quantile probability
values, $\tau_{i}$, must increase from their original quantile probability
value in order to maintain a constant multivariate quantile probability. When
the variates are independent (case 2) the individual quantile probabilities
need to be adjusted less than when the variates are maximally, negatively,
correlated using Bonferroni inequality (case 3). Therefore, the individual
quantile probability values must be adjusted slightly more for case 3 than for
case 2, since the probability of the intersection of independent events is the
product of their probabilities. Case 1 is the trivial case when the variates
are maximally, positively, correlated ($\rho=1$).

Figure 2: Top plot shows how the multivariate or joint quantile probability,
$\tau_{\mathrm{J}}$, decreases assuming a constant individual quantile
probability value, $\tau_{i}=0.90$. At ten variates, $q=10$, the Bonferroni
inequality on the lower bound on the joint quantile probability reaches zero,
$\tau_{\mathrm{J}}=0$. The Bonferroni inequality bound occurs when the
correlation among all the variates are maximally negative. However, when all
the variates are independent, the joint quantile probability only equals zero
at the limit as the number of variates reaches infinity,
$\lim\limits_{q\to\infty}(0.90)^{q}=0$. Case 1 is the trivial case when the
variates are maximally, positively, correlated ($\rho=1$). The bottom plot
shows what the individual quantile probability values must be adjusted to for
both cases in order to achieve a constant joint quantile probability value,
$\tau_{\mathrm{J}}=0.90$. The inset plot, found in the bottom plot,
demonstrates that the individual quantile probability values must be adjusted
slightly more for case 3 than for case 2.

###  2.1 Calculating the Joint Quantile Probability from a Known Distribution

Instead of bounding the quantile probability for the joint distribution, one
can numerically calculate the exact quantile probability when the variates can
be modeled from a known multivariate normal distribution. This can be done
when using concurrent univariate quantiles. The term concurrent implies that
each variate utilizes the same quantile probability, i.e.,
$\tau_{x}=\tau_{y}=\tau_{z}$. The numerical method works by utilizing
equicoordinate quantiles. To demonstrate this, an example in $\mathbb{R}^{2}$
will be used with the standard normal bivariate distribution,

$\displaystyle\begin{pmatrix}X\\\ Y\\\
\end{pmatrix}\sim\mathcal{N}\left[\boldsymbol{\mu}=\begin{pmatrix}0\\\
0\end{pmatrix},\Sigma=\begin{pmatrix}1&\rho\\\ \rho&1\end{pmatrix}\right],$
(4)

where $\rho$ denotes the correlation between the two random variables, $X$ and
$Y$. In the example provided, the marginal distributions for $X$ and $Y$ are
the just the standard normal distribution, $X\sim\mathcal{N}(0,1)$ and
$Y\sim\mathcal{N}(0,1)$. Since the parameters associated for each marginal
distribution are known, the one-sided upper tolerance bound is equivalent to
the $\tau^{\mathrm{th}}$ quantile probability, $Q(\tau=0.90)=1.2816$. Using
equations (18) and (19) while assuming an individual quantile probability
value of $0.90=\tau_{i}:i\in\left\\{X,Y\right\\}$, implies that the lower and
upper bounds for the joint quantile probability are,
$0.80\leq\tau_{{\scriptscriptstyle XY}}\leq 0.90$.

In the case where the variates are uncorrelated, $\tau_{{\scriptscriptstyle
XY}}=0.90^{2}=0.81$. The solution can be graphically seen in the second row,
middle column, of Figure 3. One can visually see that the intersection of each
univariate quantile lines and the CDF contour line occurs at the multivariate
quantile probability at a value of 0.81. For the case when $\rho=-0.99$ and
$\rho=0.99$, the intersection of each univariate quantile line and the CDF
contour lines occurs at the multivariate quantile probability at a value of
0.800 and 0.8901 respectively. The solutions are graphically shown in the
second row, first and third column, respectively, of Figure 3.

Numerical solutions can be had using the mvtnorm library in the R language
which provides a function called qmvnorm() which computes the equicoordinate
quantile function of the multivariate normal distribution for arbitrary
correlation matrices based on inversion of the probability distribution
function for the multivariate normal distribution. The algorithm utilizes a
stochastic root finding algorithm described in Bornkamp (2018). The numerical
solution can only be computed such that all the variates share the same
quantile probability, hence, equicoordinate, i.e.;
$\tau=\tau_{i}:i\in\left\\{1,\ldots,q\right\\}\ $ for $q$ number of variates.
If the above condition holds, then the numerical solution for the joint
quantile probability for arbitrary normal multivariate distributions can be
written as,

$\displaystyle g(p;\tau,C_{{\scriptscriptstyle
XY}})=\bigl{(}Q(\tau)-\text{{qmvnorm}}(p,\mathrm{corr}=C_{{\scriptscriptstyle
XY}},\mathrm{tail=``lower.tail"})\bigl{)}^{2},$ (5) $\displaystyle
G(p)=\underset{p}{\arg\min}\biggl{[}g(p;\tau,C_{{\scriptscriptstyle
XY}})\biggl{]},$ (6)

where $C_{{\scriptscriptstyle XY}}$ is the correlation matrix, tail,
corresponds to the optional input argument for the function, and $p:0\leq
x\leq 1$, is the space over all quantile probability values associated to the
multivariate CDF. Since qmvnorm() computes the equicoordinate quantile value,
this implies that the function returns just a scalar. For example, in the
bivariate case where $\rho=0$, the joint quantile probability value of $p$
that minimizes (6) is $p=0.81$. Therefore the equicoordinate quantile value
for this bivariate distribution would be the X, Y coordinate pair
$(1.2816,1.2816)$ as shown by the function output below,

$\displaystyle\
1.2816=\text{{qmvnorm}}(0.81,\mathrm{corr}=\begin{pmatrix}1&\rho\\\
\rho&1\end{pmatrix},\mathrm{tail=``lower.tail"}).$

When $\rho=0$, this is one of the few instances where one can analytically
solve for the joint quantile probability, $\tau_{{\scriptscriptstyle
XY}}=0.81$, without needing numerical methods.

Figure 3: The first row depicts the PDF contours, shown in grey, of three
standard normal bivariate normal distributions with correlation values of
$\rho\in\left\\{-0.99,0,0.99\right\\}$, one for each column. The black dash-
dotted line are the axis limits for the inset CDF contour colors below. The
second row depicts the CDF contours, shown in color, for each standard normal
bivariate distributions with different correlations. The dotted and dashed
lines correspond to the univariate ICDF, $Q(\tau_{i})$, at a quantile
probability of $\tau_{i}=0.90:i\in\left\\{X,Y\right\\}$. The diamond
represents the intersection of the two univariate ICDF lines.

###  2.2 Characterizing Uncertainty on the Joint Quantile Probability

The true covariance or correlation matrix is rarely known for data that is
multivariate normally distributed. Since there is uncertainty on the true
correlation matrix, there will also be uncertainty on the multivariate
quantile probability. We present here a methodology for characterizing the
uncertainty on the multivariate quantile probability by comparing to
concurrent univariate quantile probabilities, Algorithm 1. The methodology
utilizes a nonparametric bootstrap to estimate uncertainty by estimating the
sampling distribution on the multivariate quantile probability (Efron 1979;
Efron and Tibshirani 1986).

For each nonparametric bootstrap sample, the data is standardized.
Standardization implies that for each bootstrap sample, the bootstrap data is
centered such that its mean is the zero vector,
$\boldsymbol{\mu}=\boldsymbol{0}$. It also implies that the bootstrap data is
scaled such that the variance for each variate is unity, i.e.,
$\mathrm{diag}(S^{*})=\boldsymbol{1}_{q}$, where $S^{*}$ denotes the bootstrap
sample covariance. The standardization process implies that the bootstrap
sample correlation matrix equals the bootstrap sample covariance matrix,
$C^{*}=S^{*}$. Standardizing the bootstrap data allows one to utilize the
quantile function for the standard normal distribution for the individual
variates to calculate the random vector, $x$, such that its probability is
less than or equal to $\tau$. The standardization process alleviates the need
for characterizing uncertainty on the mean and variance, thus allowing to
exclusively characterize the affects of correlation on the multivariate
quantile probability. Algorithm 1 is thus simply resampling many possible
correlation matrices from the original data in a transformed space that allows
one to exploit the quantile function for the standard normal distribution.

In section 2 we have shown how to calculate uncertainty on the upper and lower
bound for the multivariate quantile probability when one neglects correlations
between the variates. Algorithm 1 provides a new methodology, thereby
improving the characterization of the uncertainty on the multivariate quantile
probability when the population is assumed to originate from a multivariate
distribution. If practitioners are going to use concurrent univariate
tolerance or concurrent univariate confidence intervals for the distribution’s
quantile when working with multivariate data, it is highly advised to use
algorithm 1 to assess the uncertainty on the multivariate quantile probability
and whether that level of uncertainty is acceptable. Ideally, with the methods
presented in section 4.2 and section 4.4 one should utilize the proper
multivariate statistical interval when working with multivariate data and not
revert to single variate statistical intervals.

Algorithm 1 \- R pseudo code for calculating the uncertainty on the joint
multivariate quantile probability, $\tau_{\gamma}$, by comparing to concurrent
univariate quantile probabilities, $\tau_{i}$. $\gamma$ are the
$\gamma^{\mathrm{th}}$ percentile values used to construct the simple
percentile bootstrap confidence intervals on the joint quantile probability,
$\tau_{\gamma}$. $D\in\mathbb{R}^{(n\times q)}$ is a matrix of data
originating from a $q$-dimensional multivariate distribution with $n$ samples.
Since the data is scaled with unity variance along the diagonal, this implies
the correlation matrix equals the sample covariance matrix, $C^{*}=S^{*}$.

1:$0<\tau_{i}<1,\quad 0<\gamma<1$

2:function QuantileUQ($D,\tau_{i},\gamma$)

3: $b=1000$ $\triangleright$ Bootstrap Number

4: $n=\mathrm{nrow}(D)$ $\triangleright$ Number of Samples

5: for $i$ = 1:$b$ do $\triangleright$ Bootstrap Loop

6:

7: $s_{i}=\text{{sample}}(n,\mathrm{replace=TRUE})$ $\triangleright$ Bootstrap
Indices

8: $D^{*}=D\left[s_{i},\right]$ $\triangleright$ Bootstrap Resample

9: $R^{*}=\text{{scale}}(D^{*},\mathrm{center=TRUE,scale=TRUE})$
$\triangleright$ Standardize Resample

10: $C^{*}=\text{{cor}}(R^{*})$ $\triangleright$ Bootstrap Correlation Matrix

11:
$J\left[i\right]=\underset{p}{\arg\min}\biggl{[}g(p;\tau_{i},C^{*})\biggl{]}$
$\triangleright$ (6)

12: end for

13: $\tau_{\gamma}=\text{{percentile}}(J,\gamma)$ $\triangleright$ Percentile
Bootstrap Confidence Interval

14: Return: $\tau_{\gamma}$

15:end function

##  3 Tolerance Regions in Multivariate Statistics

###  3.1 Traditional Tolerance Regions

Traditional tolerance regions can be understood as a higher dimension
generalization of the univariate two-sided tolerance intervals. Construction
of traditional tolerance regions for a $q$-variate normal distribution is
performed through Monte Carlo simulation from Algorithm 9.2 in Krishnamoorthy
and Mathew (2009) to estimate the threshold squared Mahalanobis distance, $r$.
Figure 4 depicts several examples of what traditional tolerance regions look
like for a bivariate distribution with a zero vector mean, unity variance, and
with a correlation of $\rho=0.75$. Three different proportion values, also
known as coverage values, were chosen,
$\beta\in\left\\{0.25,0.50,0.90\right\\}$ for three different sample sizes,
$n\in\left\\{10,25,1\mathrm{e}6\approx\infty\right\\}$. The tolerance regions
plotted use a constant confidence value of 90%. Larger confidence values would
increase the area of the elliptical tolerance region while smaller confidence
values would decrease the area. As the sample size increases, the elliptical
tolerance regions shrink due to increasing certainty of the multivariate
distribution’s parameters. As the sample size reaches infinity, or
equivalently when the distribution’s parameters are known, the confidence term
has no importance as there is no uncertainty on the distribution’s parameters.
There is a lack of nomenclature as to what to call these elliptical tolerance
regions when the distribution’s parameters are known; going forward, the paper
will refer to them as concentric elliptical quantiles. Elliptical tolerance
regions imply the smallest region, centered at the mean, that encompasses at
least a proportion $\beta$ of the population for a particular confidence.

One can also interpret traditional tolerance intervals through the lenses of
real analysis. Utilizing the bivariate distribution as an example, the joint
probability density function $f(x_{1},x_{2})$ is:

$\displaystyle
f(x_{1},x_{2})=\frac{\exp{\bigg{(}-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{T}\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu})}\bigg{)}}{\sqrt{(2\pi)^{2}\mathrm{det}(\Sigma)}},$
(7)

where $\mathrm{det}(\bullet)$ denotes the determinant operator and
$\boldsymbol{\mu}=[\mu_{1},\mu_{2}]$. One can write the tolerance region as,

$\displaystyle\beta=\int\int_{\Omega}f(x_{1},x_{2})\ dx_{1}\ dx_{2}:$ (8)
$\displaystyle\Omega=\left\\{\boldsymbol{x}=(x_{1},x_{2})\ |\
h(\boldsymbol{x},\bar{\boldsymbol{x}},S)\leq r\right\\},$

where $r$ is the threshold squared Mahalanobis distance. Traditional tolerance
regions imply coverage integrated outwards from the mean, $\boldsymbol{\mu}$.
This can be found directly in (21) of the form,
$(\boldsymbol{x}-\bar{\boldsymbol{x}})$. Another interpretation of (8) can be
visualized as the third column of Figure 4, when $d\leq r=4.61$, this implies
$\beta=0.90$, therefore 90% of the population resides within a squared
Mahalanobis distance of less than or equal to $r=4.61$ units. Tolerance
regions in $\mathbb{R}^{2}$ take the form of an ellipse, in $\mathbb{R}^{3}$
an ellipsoid, and for $\mathbb{R}^{q}:q>3$, the tolerance region is a hyper-
ellipsoid.

Figure 4: Various examples of normal tolerance regions for a standard normal
bivariate distribution with $\rho=0.75$. Three different proportion values are
shown, $\beta=\left\\{0.25,0.50,0.90\right\\}$ for three different sample
sizes, $n=\left\\{10,25,1\mathrm{e}6\approx\infty\right\\}$. All the tolerance
regions plotted use a constant confidence value of 90% and using algorithm 9.2
from Krishnamoorthy and Mathew (2009).

###  3.2 CDF Based Quantile Contours vs. Concentric Elliptical Quantiles for
Bivariate Normal Distributions

The easiest case for comparing concentric elliptical quantiles to CDF based
quantiles occurs for the bivariate normal distribution. Figure 5 demonstrates
how the two quantile methods compare for correlation coefficients at,
$\rho\in\left\\{-0.99,0,0.99\right\\}$, respectively for a standard normal
bivariate distribution. When $\rho=-0.99$, the concentric elliptical quantile
associated with $\beta=0.90$ are quite far from the CDF quantile contour
associated with $\tau=0.90$. An infinite number of points exists on an
ellipse. Trying to select one or a few points on the ellipse when $\rho=-0.99$
for multivariate tolerance applications is quite challenging. Maximizing a
value for one variate implies minimizing the value for another variate. For
quality control applications, testing a component at many points on a
multivariate quantile is cost prohibitive. In mechanical shock and vibration,
testing a single component at more than one point on the quantile could lead
to degradation of the component’s mechanical performance by inadvertently
fatiguing the component. When $\rho=0.99$, the concentric elliptical quantile
associated with $\beta=0.90$ extends well beyond the CDF quantile contour
associated with $\tau=0.90$ and even further past the CDF quantile contour
associated with $\tau=0.95$ (not shown). For the case where the variates are
uncorrelated, $\rho=0$, the ellipse is now a circle and is relatively close to
the CDF quantile contours. It is quite evident by comparing the CDF quantile
contours to the concentric elliptical quantiles that traditional tolerance
regions do not consider the tails of the multivariate distributions.

Figure 5: Comparing concentric elliptical quantiles to CDF based quantile
contours when the distribution parameters are known. A standard normal
bivariate distribution with various correlation coefficients at,
$\rho=\left\\{-0.99,0,0.99\right\\},$ respectively are shown.

###  3.3 CDF Based Quantile Contours vs. Simultaneous Univariate Tolerance
Intervals

When one doesn’t know the distribution’s parameters, Krishnamoorthy suggests
to use a simultaneous univariate tolerance interval to estimate the extreme
responses of a normal multivariate distribution (Krishnamoorthy and Mathew
2009). A one-sided, upper, simultaneous univariate tolerance bound, is shown
below:

$\displaystyle\tilde{T}_{\mathrm{simul},\beta}\bigg{(}1-\frac{\alpha}{q};\beta,n\bigg{)}.$
(9)

A simultaneous univariate tolerance bound is simply adjusting the confidence
parameter, $\frac{\alpha}{q}$, by a Bonferroni correction by incorporating the
number of variates, $q$, into a regular univariate tolerance interval. When
the sample size approaches infinity, the simultaneous univariate tolerance
bound converges to the standard univariate quantiles for each variate. While
extremely efficient to compute, it doesn’t provide a unique point on a
multivariate quantile but rather a point such that the coverage of the
multivariate distribution is at least $\beta$ for a particular confidence
value.

Much of this paper up until now has been designed to build up a case for how
correlation can significantly affect the CDF based quantiles for multivariate
normal distributions. All the nuances for how correlation affect the
quantiles, see Figure 15 as an example, is lost by using the Bonferroni
correction in the simultaneous univariate tolerance intervals. The
simultaneous univariate tolerance intervals do provide a desired coverage but
that particular point doesn’t satisfy the definition of a point associated to
the $\tau^{\mathrm{th}}$ quantile probability. Therefore, there is a strong
need for a new method that incorporates the correlation of the data, accounts
for the tails of the multivariate distribution, while also maintaining the
conservative properties that a point residing on the CDF based quantiles has.

##  4 A Need for a New Methodology

There are many instances where one is more concerned about extreme
multivariate responses in contrast to what squared Mahalanobis distance
encompasses a proportion of the population for a given confidence value.
Therefore, one would want a CDF based quantile compared to an elliptical
tolerance region. This is analogous to choosing a one-sided vs two-sided
statistical interval in univariate statistics.

Two main challenges exist when working with the quantile function (1) for
multivariate normal distributions. The first is that an infinite number of
vectors, $\boldsymbol{x}$, satisfy the condition $\tau\leq F(\boldsymbol{x})$.
The second, even if there was a finite set of vectors that satisfied the
condition above, no analytical solution is known for the quantile function (1)
and therefore must be numerically estimated.

The following subsections describe two methodologies for calculating
confidence intervals on multivariate normal quantiles given the challenges
listed above. The first methodology is a computationally intense Monte Carlo
meshing algorithm. The goal of the first methodology is to estimate
uncertainty on CDF based multivariate quantiles over a reasonable domain that
captures an essential, finite, subset of vectors from the quantile as they
extend towards their asymptotic limits. This methodology will be denoted
simply as confidence intervals on multivariate quantiles. Since multivariate
quantiles contain an infinite set of vectors, the goal of the second
methodology is to estimate uncertainty on the vector within the set that has
the highest probability density value.

$\displaystyle\boldsymbol{v}(\tau)=\underset{\boldsymbol{x}\ :\
\boldsymbol{x}\ \in\
Q(\tau)}{\arg\max}\biggl{[}f\bigl{(}Q(\tau)\bigl{)}\biggl{]},$ (10)

where $Q(\tau)$ is the set of vectors that are associated to the
$\tau^{\mathrm{th}}$ quantile probability and $f(\boldsymbol{x})$ is the PDF
of a $q$-variate normal distribution. The vector within the set associated to
the CDF-based quantile, $Q(\tau)$, that has the highest probability density is
given the novel term: critical point. This methodology will be denoted simply
as confidence intervals of critical points on multivariate quantiles. A
critical point can be understood as the most “likely” value on a CDF based
quantile to occur. We present without proof that a unique critical point
exists. When overlaying the CDF based quantiles with the PDF contours for
multivariate distributions that are zero mean with unity variance, the CDF
based quantiles are tangent to the PDF contours at the equicoordinate line for
all possible correlations. Therefore, any translations from the equicoordinate
line reduces the probability density.

###  4.1 Standardizing Multivariate Data

The emphasis on standardizing multivariate data or multivariate distributions
cannot be overstated for any of the methodologies presented herein. In
univariate statistics, the process of transforming data such that it has a
zero mean and unit variance is often called a z-transform, standardization, or
auto-scaling. Commonly notated as:

$\displaystyle z_{i}=\frac{x_{i}-\bar{x}}{s},$ (11)

where $z_{i}$ is the transformed observation, $x_{i}$ is the original
observation, $\bar{x}$ and $s$ are the sample mean and sample standard
deviation, respectively. Standardization for multivariate data takes an
analogous matrix form and is denoted as $\Psi(D)$,

$\displaystyle\mathcal{D}=\Psi(D):=(D-\bar{\boldsymbol{x}})\bigg{(}{\sqrt{\mathrm{diag}\big{(}{\mathrm{diag}({S})}\big{)}}}\bigg{)}^{-1}\
,$ (12)

where $\mathcal{D}$ is the transformed multivariate data and
$D\in\mathbb{R}^{(n\times q)}$ is a matrix of data originating from a
$q$-dimensional multivariate distribution with $n$ samples. The diagonal
matrix where the elements along the diagonal are the sample variances is
denotes as ${\mathrm{diag}\big{(}{\mathrm{diag}({S})}\big{)}}$ and the sample
covariance matrix is denoted as $S$. The affine transformation back to the
original domain is thus, $D=\Psi^{-1}(\mathcal{D})$. Standardization in
multivariate statistics preserves the correlation between the multivariate
data while transforming the data such that is has a zero-mean vector. The
covariance matrix of the transformed data now has unity variances along
diagonal of it, thus making the covariance matrix also the correlation matrix.

Standardization is utilized in constructing confidence intervals on
multivariate normal quantiles. When standardizing the multivariate data, only
one numerical study is required to determine the optimal mesh step size to
avoid graphical aliasing of the contours or iso-surfaces, see Section 4.2.
Constructing confidence intervals for critical points on multivariate
quantiles also requires the multivariate data to be standardized for the
numerical method to work, see Section 4.4.

###  4.2 Method 1 - Confidence Intervals on Multivariate Normal Quantiles

The purpose of this section to describe a methodology for constructing
confidence intervals on CDF based multivariate normal quantiles. The need for
this methodology occurs when the multivariate distribution’s parameters are
not known and when one requires more than just a single critical point that
resides on the quantile. Several considerations must be made to properly mesh
the support of the multivariate distribution. One, what are the bounds of the
mesh that is being tessellated. Two, what granularity is required to properly
discretize the support. To address the first point, since the multivariate
data is standardized, as is discussed in Section 4.1, we suggest that the
bounds over the domain are, $\Omega\subset\mathbb{R}^{q}=\left\\{-4\leq
x_{i}\leq 4\right\\}:i={1,\ldots,q}$ due to the low probability of values
outside of this domain. A Bonferroni lower limit for the minimum probability
over the $q$-dimensional space is calculated below for univariate
distributions up to trivariate distributions, $q=\left\\{1,2,3\right\\}$,

$\displaystyle
1-\bigg{[}q\bigg{(}1-F(4\times\boldsymbol{1}_{q})+F(-4\times\boldsymbol{1}_{q})\bigg{)}\bigg{]}=\left\\{0.999936,0.999873,0.999810\right\\},$
(13)

where $\boldsymbol{1}_{q}$ denotes a one vector of size $q$ and $F(\bullet)$
is still the CDF function for a $q$-variate normal distribution. From (13),
one can see that the domain being tessellated covers most of the codomain of
$F$.

To address the second question regarding how finely one discretizes the
support, a study was performed using three different mesh step sizes with
values of, $\left\\{0.001,0.01,0.1\right\\}$, shown in Figure 6. The study
shows little to no aliasing of the CDF based quantile contours for mesh sizes
less than 0.01 units for precision at three decimal places. Even with a mesh
size of 0.1 and a high correlation coefficient of $\rho=0.99$, the CDF based
quantile contour is nearly indistinguishable to the contours meshed at a
smaller size. When graphical aliasing does occur, it is at areas of strong
curvature, high positive correlation values, and where the PDF is highest. The
authors recommend using a mesh of 0.01 for constructing confidence intervals
for bivariate normal distributions and a mesh size of 0.1 for trivariate
normal distributions. The computational burden of meshing the support for a
trivariate normal distribution is offset using a larger mesh size of 0.1 but
at the expense of accuracy and graphical aliasing.

(a) CDF quantile contour from a bivariate distribution with a correlation of,
$\rho=-0.99$.

(b) CDF quantile contour from a bivariate distribution with a correlation of,
$\rho=0.99$.

Figure 6: CDF quantile contours represent mesh step sizes with values of,
$\left\\{0.001,0.01,0.1\right\\}$, for solid, dashed, and dotted, curves
respectively. Both plots are for a bivariate quantile probability value of
$\tau=0.90$. The solid and dashed lines are superimposed onto each other
representing that the mesh size is sufficiently small for the bivariate case.

Algorithm 2 provides the methodology to calculate confidence intervals on
trivariate normal quantiles. The algorithm can generalize to higher
dimensions, tessellation is extended over the additional dimensions.
Nevertheless, the computation burden of tessellating in higher dimensions
grows exponentially: $m^{q}$. Similarly, the algorithm can generalize to the
bivariate case and interpolation is also modified from 3 dimensions to $q=2$
dimensions. The algorithm provides two different options for bootstrapping: a
parametric and nonparametric method. Parametric bootstrapping utilizes the
original sample mean and sample covariance to randomly generate multivariate
data of the same size as the original. Efron and Hastie (2016) has suggested
that a parametric family, like a multivariate normal distribution, can act as
regularizers, smoothing out the raw data and de-emphasizing outliers. A
nonparametric bootstrap method is also provided which imposes fewer
constraints and assumptions on the underlying data. When working with the
nonparametric bootstrap approach, one must be cognizant that resampling data
with low initial sample sizes, $n$, can lead to degenerate matrices. Thus, one
must either resample again until the problem is alleviated or simply use the
parametric approach which avoids these issues.

The argument, interpolate, found in Algorithm 2, is a Boolean data type. When
the variable is true, cubic interpolation occurs over the $q$-dimensional
space, for the percentiles of the tessellated CDF mesh by an upsampling factor
of $\xi=2$. Interpolation on the tessellated mesh prior to final contour
method provides additional accuracy on the final contours that would be
originally inhibited by the memory required to store the CDF tensor at such
fine granularity. Interpolation has precedence in previous multivariate
statistical CDF work by Langrené and Warin (2021) who utilized it in their
fast multivariate empirical CDF calculations.

Algorithm 2 \- MATLAB pseudo code for a function that calculates confidence
intervals for trivariate normal CDF based quantiles. $D\in\mathbb{R}^{(n\times
3)}$ is a matrix of data originating from a $3$-dimensional multivariate
distribution with $n$ samples. $\tau$ is the desired multivariate quantile
probability value. $\gamma$ are the $\gamma^{\mathrm{th}}$ percentile values
used to construct the simple percentile bootstrap confidence intervals on the
tessellated CDF mesh denoted as, $Z_{\gamma}$. The argument, style, controls
whether to perform parametric and nonparametric bootstrapping. The argument,
interpolate, is a Boolean data type that when true, interpolates the
percentiles of the tessellated CDF mesh by an upsampling factor of $\xi$. To
denote that interpolation was performed on the tessellated CDF mesh, the breve
symbol is used, $\breve{Z}_{\gamma}$. The hat symbol on $\hat{Q}_{\gamma}$
denotes that this is an estimated subset of $Q(\tau)$, since $Q(\tau)$ is an
infinite set of vectors. Lastly, the under-bar accent in,
$\underaccent{\bar}{\hat{Q}}_{\gamma}$, denotes that the operation was
performed in the standardized domain. The inverse, affine, transformation of
the quantile back to the original domain of the data is simply denoted without
the under-bar.

1:$0<\tau<1,\quad 0<\gamma<1,\quad$

2:function TrivariateQuantile($D,\tau,\gamma$, style, interpolate)

3: $b=1000$ $\triangleright$ Bootstrap Number

4: meshStep = 0.1

5: $[n,q]=\text{{size}}(D)$ $\triangleright$ Number of Samples and Variates

6: $\mathcal{D}=\text{{normalize}}(D,\mathrm{``zscore"})$ $\triangleright$
Standardize Data

7: $S=\text{{cov}}(\mathcal{D})$ $\triangleright$ Covariance/Correlation
Matrix

8: x = [-4 : meshStep : 4]

9: y = [-4 : meshStep : 4]

10: z = [-4 : meshStep : 4]

11: [X,Y,Z] = meshgrid(x,y,z)

12: XYZ = [X(:), Y(:), Z(:)]

13: MC = nan(length(y),length(x),length(z),b) $\triangleright$ Pre-Allocate
nan tensor

14: for $i$ = 1:$b$ do $\triangleright$ Monte Carlo Loop

15: if style == ‘parametric’ then $\triangleright$ Parametric Bootstrap

16:
$\mathcal{D}^{*}=\text{{mvnrnd}}(n,\mathrm{mu}=\boldsymbol{0},\mathrm{Sigma}=S)$
$\triangleright$ Parametric Resample

17: else if style == ‘nonparametric’ then $\triangleright$ Nonparametric
Bootstrap

18: $\mathcal{D}^{*}=\text{{datasample}}(\mathcal{D},n,\mathrm{replace=TRUE})$
$\triangleright$ Nonparametric Resample

19: end if

20: $\bar{\boldsymbol{x}}^{*}=\text{{mean}}(\mathcal{D}^{*},\mathrm{DIM}=1)$
$\triangleright$ Bootstrap sample Mean

21: $S^{*}=\text{{cov}}(\mathcal{D}^{*})$ $\triangleright$ Bootstrap sample
covariance

22: $Z^{*}=\text{{mvncdf}}(\mathrm{XYZ},\bar{\boldsymbol{x}}^{*},S^{*})$
$\triangleright$ Bootstrap CDF

23: $\mathrm{MC}(:,:,:,i)=Z^{*}$; $\triangleright$ Store Bootstrap CDF

24: end for

25:
$Z_{\gamma}=\text{{percentile}}(\mathrm{MC},\gamma,\mathrm{DIM}=q+1)\textbf{}$
$\triangleright$ percentile the CDF Mesh

26: if interpolate == TRUE then

27:
$[\breve{Z}_{\gamma},\mathrm{x},\mathrm{y},\mathrm{z}]=\text{{interp3}}(Z_{\gamma},\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{\xi}=2)$

28:
$\underaccent{\bar}{\hat{Q}}_{\gamma}=\text{{contourc}}(\mathrm{x},\mathrm{y},\mathrm{z},\breve{Z}_{\gamma},\tau)$
$\triangleright$ Extract Contours

29: else

30:
$\underaccent{\bar}{\hat{Q}}_{\gamma}=\text{{contourc}}(\mathrm{x},\mathrm{y},\mathrm{z},Z_{\gamma},\tau)$
$\triangleright$ Extract Contours

31: end if

32: ${\hat{Q}}_{\gamma}=\Psi^{-1}(\underaccent{\bar}{\hat{Q}}_{\gamma})$
$\triangleright$ Inverse Transform to Original Domain

33: Return: ${\hat{Q}}_{\gamma}$ $\triangleright$ Percentile Bootstrap
Confidence Interval

34:end function

(a) Pictograph describing how to generate bivariate normal CDF based quantile
contours. Each bootstrap within the Monte Carlo simulation will create a
tessellated grid with dimensions $m\times m$ for which the CDF is calculated
over. From Algorithm 2, the support is tessellated mesh -4 to 4 with a mesh
step size of 0.01, thus $m=801$. After all $b$ Monte Carlo simulations have
been performed, the CDF grids are now a tensor with dimensions $m\times
m\times b$. The $\gamma^{\mathrm{th}}$ percentile is performed along the third
dimension of the tensor, resulting back to a final dimension $m\times m$. The
$\gamma^{th}$ percentile provides uncertainty on the CDF values from all the
various multivariate bootstrapped data. It is important to note that the
$\gamma^{\mathrm{th}}$ percentile is over the tessellated CDF tensor, which is
not to be confused with the $\tau^{\mathrm{th}}$ quantile probability which
refers to the quantile probability associated with the multivariate normal
distribution.

(b) Once the $\gamma^{\mathrm{th}}$ percentile is computed over the
tessellated CDF tensor from 7(a), the contours associated to the
$\tau^{\mathrm{th}}$ quantile probability can be estimated. The pictograph
above provides a rudimentary example for a multivariate normal quantiles
associated with $\tau=0.90$ highlighted in black with white colored font. The
uncertainty associated with the set of vectors comprising a multivariate CDF
quantile is denoted as ${\hat{Q}}_{\gamma}(\tau)$. The blank spots between the
bottom left and top right of the grid would represent monotonically increasing
CDF values. A similar procedure can be performed for trivariate normal
distributions, however, the quantile associated for ${\hat{Q}}_{\gamma}(\tau)$
is an iso-surface from a tessellated cube that with dimensions $m\times
m\times m$.

Figure 7: Pictographs describing the process of generating confidence
intervals for CDF based bivariate normal quantiles,
${\hat{Q}}_{\gamma}(\tau)$. The methodology is similar for higher dimensions
such as trivariate normal quantiles, but with much higher computational burden
for an increase in dimensionality.

In Figure 8 we compare the bivariate version of algorithm 2 to the only
current methodology for generating multivariate quantiles, the latter is only
applicable in the bivariate case (Klein and Kneib 2020). Fifty random
bivariate vectors were sampled from the same bivariate distribution to create
a synthetic data set, $D\in\mathbb{R}^{(50\times 2)}$, shown as blue circles.
The $i^{\mathrm{th}}$ random vector are sampled from a bivariate normal
distribution with parameters,
$D_{i}\sim\mathcal{N}\bigg{(}\boldsymbol{\mu}=[6,10]),\Sigma=\begin{bmatrix}1&0.50\\\
0.50&1\end{bmatrix}\bigg{)}$. The true contour, $\tau=0.90$, is shown as a
thick black line. The directional bivariate quantile from Klein and Kneib
(2020) is shown as a red dotted line. The methodology described above for
constructing confidence intervals on multivariate normal quantiles is shown as
the shaded grey region for the quantile probability of $\tau=0.90$ using a
two-sided confidence interval value of 95%. To generate the quantiles,
$b=5000$, was used as the number of Monte Carlo simulations followed by using
the bias corrected and accelerated method (BCa) bootstrap confidence interval
method. The BCa bootstrap confidence interval method was utilized instead of
the basic percentile method because it generally has smaller coverage error
(Davison and Hinkley 1997; Carpenter and Bithell 2000). The algorithms
presented here use the basic percentile method because implementing the (BCa)
methods are more complicated. The upper confidence interval is shown as a line
with diamond markers, while the lower confidence interval is shown as a line
with square markers. Figure 8 demonstrates that even with 50 samples, the
current directional bivariate method provides a rough empirical approximation
to the true bivariate contour but with no measure of uncertainty. The
methodology described here provides not only uncertainty on the distribution’s
quantile but also the smoothness one expects. All the quantiles shown in
Figure 8 asymptotically extend towards infinity with respect to each variate,
yet one can see that domain for which the support is meshed captures the
essence of the underlying bivariate quantiles.

Figure 8: Bivariate normal distribution example comparing the directional
bivariate quantile methodology and the present methodology described above for
constructing confidence intervals on multivariate normal quantiles. The true
contour, $Q(\tau=0.90)$, is shown as a thick black line. The directional
bivariate quantile from Klein and Kneib (2020) is shown as a red dotted line.
The methodology described above for constructing confidence intervals on
multivariate normal quantiles, ${\hat{Q}}_{\gamma}(\tau)$, are shown as the
shaded grey region for the quantile probability of $\tau=0.90$ using a two-
sided confidence interval of 95%. To generate the quantiles, $b=5000$, was
used as the number of Monte Carlo simulations followed by using the BCa
bootstrap confidence interval method. The upper confidence interval,
${\hat{Q}}_{0.975}(\tau)$, is shown as a line with diamond markers, while the
lower confidence interval, ${\hat{Q}}_{0.025}(\tau)$, is shown as a line with
square markers.

The computational time on a 12th Gen Intel® Core™ i9-12950HX 2.30 GHz
processor is about 189 seconds when using 5000 bootstrap resamples for the
present methodology compared to about 0.90 seconds for the directional
bivariate quantile method. The smoothness of the calculated quantile, the
extended domain for which the quantile covers, and the measure of uncertainty
on the quantile does come at a computational cost of approximately two orders
of magnitude, nonetheless, the overall computational time is still reasonable
for most applications. Computation time for iso-surfaces, representing
quantiles for trivariate distributions, require significantly more computation
time of approximately 6 hours for 1000 bootstrap resamples. As of to-date,
algorithm 2 is the only method for calculating confidence intervals for
quantiles associated with trivariate normal distributions.

###  4.3 Critical Point Selection

The single point on a multivariate quantile that has the maximum relative
likelihood, or equivalently, highest probability density of occurring, is
known as the critical point. The critical point will be denoted as the bold
vector, $\boldsymbol{v}(\tau)$. These points can be estimated from the finite
subset of vectors comprising the multivariate quantiles using a variety of
techniques. If the data has been standardized, the critical point can be
represented as $\boldsymbol{v}=\boldsymbol{1}_{q}v$ where $v$ is the
equicoordinate value of the CDF based quantile, i.e., $v=x_{1}=\ldots=x_{q}$.
In $\mathbb{R}^{2}$ the intersection of the equation $y=x$ with the CDF based
quantile, represented as a contour, is the critical point. In $\mathbb{R}^{3}$
the intersection of the parametric equation:

$\displaystyle\begin{cases}x(t)=t\\\ y(t)=t\\\ z(t)=t\end{cases},$ (14)

with the CDF based quantile, represented now as an iso-surface, is the
critical point.

Figure 9 demonstrates the critical points, $\boldsymbol{v}(\tau)$, from a
standardized bivariate normal distribution for three different quantile
probabilities, $\tau=\left\\{0.01,0.5,0.99\right\\}$, and three different
correlation values, $\rho=\left\\{-0.99,0,0.99\right\\}$. Since the
distribution parameters are known for this example, one can select the vector
that is associated with the largest probability density which is simply the
vector from the quantile that has the maximum PDF value. Similarly, one can
translate the desired CDF quantile from data that has been standardized into
the first quadrant, followed by finding the vector that minimizes the $L^{2}$
norm. Translating the CDF quantiles into the first quadrant, or in higher
dimensions orthant, forces all of the elements within the set of vectors that
comprises a quantile to be positive. One can also see that the vector on the
multivariate quantile that minimizes the squared Mahalanobis distance is
generally the critical point but not always as for the case when $\tau=0.01$
and $\rho=0.99$, shown as a the purple diamond marker in the top right plot in
Figure 9. Regardless of the method used to calculate the critical point, if
the data is standardized, any of the methods, besides the minimum squared
Mahalanobis distance, yield the same answer. When the parameters of the
distribution are not known, either the translation or intersection method can
be used when working with a particular quantile. To denote uncertainty on the
critical points, $\boldsymbol{v}_{\gamma}(\tau)$, denotes the
$\gamma^{\mathrm{th}}$ percentile on a critical point from any of the
bootstrap confidence intervals associated for the $\tau^{\mathrm{th}}$
multivariate quantile probability.

Figure 9: Critical points, $\boldsymbol{v}(\tau)$, from a standardized
bivariate normal distribution for three different quantile probabilities,
$\tau=\left\\{0.01,0.5,0.99\right\\}$, shown along the $x_{2}$ axis, and three
different correlation values, $\rho=\left\\{-0.99,0,0.99\right\\}$, shown
along the $x_{1}$ axis.

###  4.4 Method 2 - Confidence Intervals of Critical Points on Multivariate
Quantiles

Under certain circumstances, one may not need a finite, subset of vectors for
the multivariate quantiles but rather just the critical point for each
quantile. When this is the case, one can substantially speed up the
computation time by directly calculating only the critical point for a
particular quantile, which is the topic for this section. One can utilize the
mvtnorm library in the R language which provides a function called qmvnorm to
help compute the critical values for multivariate normal distribution.

When the parameters of the multivariate normal distribution are known, thus
confidence is irrelevant, one must first calculate the correlation matrix,
$C$, from the known covariance matrix, $C=R\Sigma
R:R=\bigg{(}\sqrt{\text{diag}(S)}\bigg{)}^{-1}$. One then utilizes qmvnorm
with the known correlation matrix to solve for the equicoordinate value, $v$,
for a standard multivariate normal distribution. When the multivariate normal
distribution is standardized, zero mean and unity variance, the equicoordinate
vector can then be transformed to any arbitrary domain for all multivariate
normal distributions with known distribution parameters,
${\boldsymbol{v}}=\Psi(v)=\bigg{(}\boldsymbol{1}_{q}v\times\sqrt{\text{diag}(\Sigma)}\bigg{)}+\boldsymbol{\mu}$
.

When the parameters of the multivariate normal distribution are not known, a
Monte Carlo bootstrap procedure can be performed on the data. Algorithm 3
provides an example of how to construct confidence intervals on the critical
points on multivariate quantiles. Regardless of whether the parameters of the
multivariate normal distribution are known or not, the calculations for the
equicoordinate vector should be performed on a standard multivariate normal
distribution. The critical point is then constructed by transforming the
equicoordinate value to the domain of any arbitrary multivariate normal
distribution.

Algorithm 3 \- R pseudo code for calculating confidence intervals on the
critical points on multivariate quantiles,
$\hat{\boldsymbol{v}}_{\gamma}(\tau)$. $\gamma$ are the $\gamma^{\mathrm{th}}$
percentile values used to construct the simple percentile bootstrap confidence
intervals on the critical points for a particular quantile probability,
$\tau$. $D\in\mathbb{R}^{(n\times q)}$ is a matrix of data originating from a
$q$-dimensional multivariate distribution with $n$ samples.

1:$0<\tau<1,\quad 0<\gamma<1$

2:function CriticalPointUQ($D,\tau,\gamma$)

3: $b=1000$ $\triangleright$ Bootstrap Number

4: $n=\mathrm{nrow}(D)$ $\triangleright$ Number of Samples

5: for $i$ = 1:$b$ do $\triangleright$ Bootstrap Loop

6:

7: $s_{i}=\text{{sample}}(n,\mathrm{replace=TRUE})$ $\triangleright$ Bootstrap
Indices

8: $D^{*}=D\left[s_{i},\right]$ $\triangleright$ Bootstrap Resample

9: $\mathcal{D}^{*}=\text{{scale}}(D^{*},\mathrm{center=TRUE,scale=TRUE})$
$\triangleright$ Standardize Resample

10: $C^{*}=\text{{cor}}(\mathcal{D}^{*})$ $\triangleright$ Bootstrap
Correlation Matrix

11:

12:
$v^{*}=\text{{qmvnorm}}(\tau,\mathrm{corr=}C^{*},\mathrm{tail=``lower.tail"})$

13:

14: $\mu^{*}=\text{{apply}}(D^{*},2,\text{{mean}})$ $\triangleright$ Bootstrap
mean

15: $\Sigma^{*}=\text{{apply}}(D^{*},2,\text{{sd}})$ $\triangleright$
Bootstrap covariance

16: $q=\text{{dim}}(D^{*})[2]$ $\triangleright$ Number of Variates

17: $V^{*}[i,]=\big{(}\text{{rep}}(v^{*},q)\times\Sigma^{*}\big{)}+\mu^{*}$
$\triangleright$ Critical Point in Original Domain

18: end for

19: $\hat{\boldsymbol{v}}_{\gamma}(\tau)=\text{{percentile}}(V^{*},\gamma)$
$\triangleright$ Percentile Bootstrap Confidence Interval

20: Return: $\hat{\boldsymbol{v}}_{\gamma}(\tau)$

21:end function

##  5 Results

###  5.1 Numerical Simulation

####  5.1.1 Algorithm 1 - Characterizing Uncertainty on the Multivariate
Quantile Probability

A numerical simulation was performed to show the efficacy of Algorithm 1.
Recall that Algorithm 1 characterizes the uncertainty on the multivariate
quantile probability, $\tau$, by comparing to concurrent univariate quantile
probabilities, $\tau_{1}=\ldots=\tau_{q}$. Two variables were investigated,
the number of variates, $q$, and, $n$, the sample size of the data. The number
of variates ranged from 2 to 4 (inclusively), $2\leq q\leq 4$, and three
different sample sizes were utilized in the simulation,
$n\in\left\\{30,300,3000\right\\}$, thus leading to 9 combinations to test. A
hundred Monte Carlo simulations were performed for each combination of
variables listed in Table 1. For each Monte Carlo simulation, a random
correlation matrix, $C_{q}$,was generated using the $C$-vine method from
Lewandowski et al. (2009). The $C$-vine method requires sampling from a beta
distribution which controls how much the distribution of partial correlations
is concentrated around ±1 using just one parameter, $c$. For the context of
the simulation, $c=2$, which provided correlation matrices with fairly large
off-diagonal values. A random mean vector of dimension $q$ was sampled from a
continuous uniform distribution: $\boldsymbol{\mu}\sim\boldsymbol{1}_{q}\
\mathcal{U}(0,100)$, where $\mathcal{U}(a,b)$ denotes a continuous uniform
distribution over the support $\left[a,b\right]$. For each Monte Carlo
simulation, data was randomly sampled from a multivariate normal distribution,
$D\sim\mathrm{mvnrnd}(\boldsymbol{\mu},C_{q},n)$, where $D$ denotes a matrix
of $n$ random vectors chosen from the same multivariate normal distribution,
with mean vector $\boldsymbol{\mu}$ and correlation matrix $C_{q}$.

Within each of the one hundred Monte Carlo simulations, the true joint
coverage, $\tau_{\mathrm{true}}$, was calculated using (6). Similarly, $D$, is
passed into Algorithm 1 to calculate a two-sided uncertainty interval for the
multivariate quantile probability with a desired confidence of $95\%$. The
quantile probability for each variate was chosen to be
$\tau_{i}=0.90:i\in\left\\{1,\ldots,q\right\\}$. Within Algorithm 1, the same
1000 bootstrap resamples were utilized, however, the bootstrap confidence
interval method used was altered from the basic percentile to the BCa
percentile method. The 95% confidence interval on $\tau_{\mathrm{true}}$ is
denoted as $\left[\tau_{{\scriptscriptstyle 0.025}},\tau_{{\scriptscriptstyle
0.975}}\right]$. Finally, at the end of each Monte Carlo simulation, if the
true multivariate coverage was bounded by the uncertainty calculated by
Algorithm 1, $\tau_{{\scriptscriptstyle
0.025}}\leq\tau_{\mathrm{true}}\leq\tau_{{\scriptscriptstyle 0.975}}$, a value
of true was returned, else false. The final column of Table 1 tabulates the
calculated confidences for each combination of variables performed in the
simulation. Overall, for numerous sample sizes and variates, Algorithm 1
provides an overall confidence of $94\%$ which is very close to the true
expected confidence value of $95\%$.

$q$ | $n$ | $p$
---|---|---
2 | 30 | 0.93
2 | 300 | 0.95
2 | 3000 | 0.89
3 | 30 | 0.94
3 | 300 | 0.98
3 | 3000 | 0.93
4 | 30 | 0.92
4 | 300 | 0.95
4 | 3000 | 0.93
AVE | 0.94
Table 1: Numerical simulation results from 100 Monte Carlo simulations for
each combination of variables on the joint multivariate quantile probability
with a desired confidence of $95\%$. Each individual variate had a quantile
probability with the associated value of,
$\tau_{i}=0.90:i\in\left\\{1,\ldots,q\right\\}$.

####  5.1.2 Algorithm 2 - Confidence Intervals on Multivariate Normal
Quantiles

Designing a numerical simulation to test Algorithm 2 is not a trivial task due
to the dimensionality and complexity of the problem space. At this current
moment, we believe that the current computation burden prevents running a
simulation study directly for confidence intervals on multivariate quantiles
for three variates. In lieu of this, the simulation study will use a modified
version of Algorithm 2 designed for two variates. There are numerous
combinations of variables that must be tested against to verify that that the
methodologies provide the correct solution for a particular confidence value.
For bivariate distributions these variables are over the space of all possible
correlation coefficients, $-1\leq\rho\leq 1$, quantile probabilities,
$0\leq\tau\leq 1$, and a host of possible sample sizes. Normalization of the
data allows one to avoid testing over the space of all possible mean vectors
and variances. Figure 10 demonstrates the complexity in verifying confidence
intervals on CDF based multivariate quantiles. The example given is for a
bivariate distribution with arbitrary parameters shown in the plot title. The
red dashed contour is the ground truth, $Q(\tau)$, for this bivariate
distribution with $\tau=0.4185$. The solid blue contour,
$\hat{Q}_{0.025}(\tau)$, and the dotted yellow contour,
$\hat{Q}_{0.975}(\tau)$, denote the lower and upper confidence intervals for a
two-sided 95% confidence interval for the multivariate quantile for an
arbitrary sample size of $n=28$. One can see that the true quantile contour
temporarily intersects the blue contour representing the lower confidence
bound on the quantile contour. Ideally, the lower and upper confidence bounds
on the $\tau^{\mathrm{th}}$ quantile would “sandwich” the ground truth
contour. When running Monte Carlo simulation studies for confidence intervals
on multivariate quantiles, one could classify any departure from the
confidence bounds as false. However, this classification is exceedingly
conservative as it foregoes most of the contour. Since the contours are
technically an infinite set of points that satisfy (1), any notion of
proportion or “majority” in this context is an ill-posed statement. Therefore,
running a simulation in this manner would provide erroneous $p$-values for the
presented method.

Figure 10: The red dashed contour is the ground truth, $Q(\tau)$, for this
bivariate distribution for a quantile probability of $\tau=0.4185$. The solid
blue contour, $\hat{Q}_{0.025}(\tau)$, and the dotted yellow contour,
$\hat{Q}_{0.975}(\tau)$, denote the lower and upper confidence intervals for
the multivariate quantile for an arbitrary sample size of $n=28$. One can see
that the true quantile temporarily intersects the blue contour representing
the lower confidence bound on the quantile.

The authors propose a simulation based on the proportion of the multivariate
population that the multivariate quantile bounds; the bivariate case is shown
below:

$\displaystyle\beta=\int\int_{\Omega}f(x_{1},x_{2})\ dx_{1}\ dx_{2}:$
$\displaystyle\Omega=\left\\{(x_{1},x_{2})\ |\ \boldsymbol{x}\leq
Q(\tau)\right\\}.$ (15)

A bivariate numerical simulation study was performed utilizing $\hat{\beta}$.
Three different correlation coefficients, $\rho=\left\\{-0.9,0,0.9\right\\}$,
and quantile probabilities, $\tau=\left\\{0.7,0.8,0.9\right\\}$, were used.
For each combination of the correlation coefficients and quantile
probabilities, $\hat{\beta}$, was computed as the ground truth using
$Q(\tau)$. The numerical simulation study utilizes 2 different sample sizes,
$n=\left\\{30,300\right\\}$ to test the efficacy of the methodology for
constructing confidence intervals on CDF based multivariate quantiles. For
each combination of $n$, 250 Monte Carlo simulations were performed where a
sample size of $n$ was drawn from the true bivariate distribution. A one-
sided, upper, confidence interval on the multivariate quantile was performed
for that particular sample at a confidence value of 95%. The upper bound of
the confidence interval was utilized to calculate
$\hat{\beta}_{i}:i=1,\ldots,250$ using the same one hundred thousand random
multivariate samples. Finally, after all 250 Monte Carlo simulations were
performed, a $p$-value was estimated,
$\mathrm{mean}(\hat{\beta}_{i}\geq\hat{\beta})$. This procedure was carried
out for all the combinations of variables, shown as the three columns, in
Table 2. The number of Monte Carlo simulations may seem small with 250, but
within each of those simulations, 1000 Monte Carlo bootstraps are performed to
construct the contours for the outer Monte Carlo simulation. Therefore, the
simulation study presented is very computationally time consuming, taking 2-3
days to compute. Hence, the authors chose to test the bivariate case of
Algorithm 2 which differs only in the tessellation of the support for the
distribution.

One can see that the $p$-values calculated are highly correlated with the
correlation coefficient, shown as the last column of Table 2. When $\rho=-0.9$
the $p$-value is much greater than the theoretical value of $p=0.95$, however,
the method is maintaining a conservative estimate. When $\rho\geq 0$, one can
see that the $p$-values are close to the expected theoretical value of
$p=0.95$. The methodology does not appear to have a bias towards the smaller
sample size of $n=30$ compared to when $n=300$ due to using the BCa bootstrap
confidence interval method. The results provided here confirm the efficacy
methodology at least in the bivariate case. For higher dimension cases, if the
support can be tessellated sufficiently, then the method should extrapolate to
accommodate them.

$\rho$ | $\tau$ | $n$ | $p$ | average
---|---|---|---|---
-0.9 | 0.7 | 30 | 0.984 | 0.997
300 | 1
0.8 | 30 | 1
300 | 1
0.9 | 30 | 1
300 | 1
0 | 0.7 | 30 | 0.960 | 0.951
300 | 0.956
0.8 | 30 | 0.940
300 | 0.956
0.9 | 30 | 0.932
300 | 0.964
0.9 | 0.7 | 30 | 0.944 | 0.944
300 | 0.941
0.8 | 30 | 0.932
300 | 0.964
0.9 | 30 | 0.948
300 | 0.936
average | 0.964
Table 2: A bivariate numerical simulation study was performed utilizing
$\hat{\beta}$. Three different correlation coefficients,
$\rho=\left\\{-0.9,0,0.9\right\\}$, and quantile probabilities,
$\tau=\left\\{0.7,0.8,0.9\right\\}$, were used. For each combination of the
correlation coefficients and quantile probabilities, $\hat{\beta}$, was
computed as the ground truth using $Q(\tau)$. The numerical simulation study
utilizes 2 different sample sizes, $n=\left\\{30,300\right\\}$ to test the
efficacy of the methodology for constructing confidence intervals on CDF based
bivariate quantiles.

####  5.1.3 Algorithm 3 - Confidence Intervals of Critical Points on
Multivariate Quantiles

Designing a simulation study for Algorithm 3 requires a different approach due
to the nature of the problem space. There is ambiguity with how to classify
whether the true critical point is within the uncertainty interval for the
critical points. Ideally,
$\hat{\boldsymbol{v}}_{\frac{\alpha}{2}}(\tau)\leq\boldsymbol{v}(\tau)\leq\hat{\boldsymbol{v}}_{1-\frac{\alpha}{2}}(\tau)$,
but one must pick a metric for measuring this criterion because the critical
points are not scalars but vectors. For example, if one element in the
uncertainty interval for the critical point vector doesn’t satisfy this
criteria, it may be too conservative to reject the entire critical point. To
avoid having to pick an arbitrary metric, one can utilize the fact that the
CDF of the true critical point provides the $\tau^{\mathrm{th}}$ quantile
probability, $F(\boldsymbol{v}(\tau))=\tau$. The corollary is that one can now
assess the probability whether the CDF value of the one-sided, $p=95\%$,
confidence interval on the critical point is at least greater or equal to the
$\tau^{\mathrm{th}}$ quantile probability. One can therefore verify that the
confidence interval on the critical point is applicable to the correct
quantile probability. It is assumed that the confidence interval on the
critical point still maintains the property that it is the point associated
with the $\tau^{\mathrm{th}}$ quantile probability with highest PDF value.
With very large sample sizes, the confidence interval on the critical point
does converge to the theoretical critical point.

Therefore a numerical simulation study was performed for three different
number of variates, $q\in\left\\{2,3,4\right\\}$, two different sample sizes,
$n\in\left\\{30,300\right\\}$, and three different quantile probabilities,
$\tau\in\left\\{0.7,0.8,0.9\right\\}$. For each combination of variables, 250
Monte Carlo simulations were performed. For each Monte Carlo simulation, $n$
number of random vectors were drawn from a standardized $q$-variate normal
distribution with a random correlation matrix, $C_{q}$, using the $C$-vine
method where $c=2$. For each Monte Carlo simulation, a one-sided, upper,
confidence interval on the critical point for the $\tau^{\mathrm{th}}$
multivariate quantile probability was calculated, denoted as,
$\boldsymbol{v}_{(1-\alpha,i)}:i=1,\ldots,250$, where $i$ denotes the
$i^{\mathrm{th}}$ Monte Carlo simulation. The CDF value of the the critical
point thus forms a one-sided, upper, confidence bound on the
$\tau^{\mathrm{th}}$ quantile probability, denoted as,

$\displaystyle\tau_{(1-\alpha,i)}=F\bigg{(}\hat{\boldsymbol{v}}_{(1-\alpha,i)}\bigg{)}.$
(16)

Finally, after all 250 Monte Carlo simulations were performed, a $p$-value was
estimated using (16) as,

$\displaystyle
P\Big{(}\tau_{(1-\alpha,i)}\geq\tau\Big{)}=\mathrm{mean}\Big{(}\tau_{(1-\alpha,i)}\geq\tau\Big{)}.$
(17)

This procedure was carried out for all the combinations of variables, shown as
the three columns, in Table 3. For a one-sided confidence interval with 95%
confidence one expects the $p$-values to be close to 0.95 but the simulation
results were 0.984. The default convergence for options for the qmvnorm() is a
maximum of 500 iterations or tolerance on the quantile of 0.001. The
simulation results provided here increased this default value to a maximum of
5000 iterations which improved the average $p$-value for the entire simulation
from 0.993 to 0.984. The results from the numerical simulation in Table 3
provides robust evidence that the numerical methods provided in Section 4.4
for constructing confidence intervals on critical points residing on
multivariate quantiles are slightly more conservative than their theoretical
confidence value of $p=0.95$.

$q$ | $n$ | $\tau$ | $p$
---|---|---|---
2 | 30 | 0.7 | 0.980
0.8 | 0.964
0.9 | 0.948
300 | 0.7 | 0.980
0.8 | 0.996
0.9 | 0.996
3 | 30 | 0.7 | 0.976
0.8 | 0.964
0.9 | 0.968
300 | 0.7 | 0.988
0.8 | 0.988
0.9 | 0.988
4 | 30 | 0.7 | 0.992
0.8 | 0.992
0.9 | 0.988
300 | 0.7 | 1
0.8 | 1
0.9 | 1
average | 0.984
Table 3: A numerical simulation study was performed for three different
variates, $q\in\left\\{2,3,4\right\\}$, two different sample sizes,
$n\in\left\\{30,300\right\\}$, and three different quantile probabilities,
$\tau\in\left\\{0.7,0.8,0.9\right\\}$. 250 Monte Carlo simulations were
performed generating a $p$-value for each combination of variables. A
confidence value of 95% was used throughout the study.

###  5.2 Shock Data Case Study for Environmental Specification

Over the past decade, there has been strong research interest in multi-axis
testing as opposed to single-axis testing within the field of mechanical shock
and vibration. Yet, many test standards still currently utilize only single-
axis shock and vibration testing such as NAVMAT P-9492 from the U.S. Navy or
the electronics industry’s JESD22-B103B (noa a, b, c). The benefits of multi-
axis testing are that it excites all modes simultaneously and induces a more
realistic stress loading condition while alleviating the issue where the order
in which orthogonal uni-axial excitation is applied has a significant effect
on fatigue failure (Whiteman and Berman 2002). Sisemore et al. (2020) provides
information about the development of multi-axis shock test specifications
while brashly stating, “writing a three-axis shock test specification is not
particularly complex”. This statement may be true when one knows the true
underlying shock excitation or its respective shock response spectra. In
practice, one rarely knows the exact environmental vibration or shock
excitation levels. Nelson (2017) provides four different derivation methods to
develop six degree of freedom (6-DOF) input specifications for multi-input,
multi-output (MIMO) mechanical vibration testing. The current literature
neglects to mention about statistical conservatism and how it applies to
multi-axis shock and vibration testing.

Environmental specification implies deriving what are the extreme shock or
vibration environmental levels that a system or component may experience in
the field in a statistical sense given some specified risk tolerance. For
single-axis shock and vibration testing this amounts to using the collected
time series acceleration field data from tri-axial accelerometers and
transforming the data into the spectral domain; power spectral density (PSD)
or shock response spectrum (SRS) curves for shock and vibration data
respectively (Lalanne 2014). At each frequency value in the spectral domain,
the data are assumed to be either normally or log-normally distributed. The
Shapiro-Wilk or Anderson Darling test can be used to test for normality. Most
engineers who are qualifying a system or component in the domain of shock and
vibration are concerned about the extreme responses across a band/range of
frequencies, thus a one-sided, upper, normal tolerance bound is often
constructed to bound at least a proportion of the data for a particular
confidence value. The case study provide here will demonstrate how to derive a
three-axis environmental specification for mechanical shock from field data
that is statistically conservative. A fully defined 6-DOF vibration test is a
$6\times 6$ spectral density matrix, the diagonal includes the three
translation directions X, Y, and Z and the rotations associated for each
spatial direction (Jacobs-Omalley 2014; Nelson and Jacobs-Omalley 2015). For
simplicity, the environmental specification here will be designed for 3-DOF
shock test for just the translation directions.

Environmental field data, time series acceleration, was collected from a test
body instrumented with a tri-axial accelerometer. The test body was exposed to
the same environmental shock nine times, thus creating an ensemble of tri-
axial accelerometer data. Each time series was then transformed into the
spectral domain forming the shock data set, SRS in the case for mechanical
shock. The SRS calculation utilized the maxi-max absolute acceleration (MMAA)
method with a 5% damping factor (Lalanne 2014). The first column of Figure 11
plots the SRS for an ensemble with each curve representing a SRS from each
sample. The rows indicate the SRS ensemble for each spatial direction. At each
frequency in the SRS domain, one would compute a pointwise, three-axis,
environmental specification for mechanical shock. To illustrate this pointwise
procedure, the SRS data occurring at an arbitrary frequency of 200.24 Hz
(shown as a dashed line for all the plots in the first column) will be further
analyzed. The second column of Figure 11 are swarm plots representing the SRS
values at a frequency of 200.24 Hz; each row still corresponds the SRS values
in the X, Y, or Z spatial directions respectively.

Figure 11: Tri-axial shock data in the shock response spectrum (SRS) domain using the maxi-max absolute acceleration (MMAA) method with a 5% damping factor. Each row corresponds the SRS MMAA for each spatial direction X, Y, and Z respectively. The first column represents the SRS over the entire frequency domain of interest from 20 Hz to 2000 Hz. The second column are swarm plots of the SRS values at a frequency of 200.24 Hz with random horizontal jitter to avoid overlapping data points for visualization purposes. The nine SRS values for each channel direction will comprise the columns of the data set, $D\in\mathbb{R}^{(9\mathrm{x}3)}$. X | Y | Z | $d_{j}$
---|---|---|---
8.49 | 5.76 | 2.75 | 4.99
6.44 | 7.81 | 3.80 | 2.12
5.26 | 5.65 | 2.67 | 1.25
3.27 | 4.27 | 3.27 | 3.36
4.81 | 5.65 | 2.47 | 2.33
3.66 | 10.64 | 3.24 | 1.95
4.96 | 6.63 | 3.62 | 1.07
5.23 | 14.61 | 3.39 | 4.64
5.27 | 10.14 | 4.06 | 2.28
Table 4:  The first three columns of the table are the tabulated data from the
second column of the swarm plots in Figure 11 for the three spatial directions
X, Y and Z respectively. Each row of the first three columns is considered as
one multivariate observation. The final column of the table is the squared
Mahalanobis distance, $d_{j}$, for each multivariate observation and is
calculated using (21) for the purpose of checking the multivariate normality
assumption.

For convenience, the SRS values occurring at 200.24 Hz for each spatial
direction are also tabulated as the first three columns in Table 4. Nine SRS
values for each of the three spatial directions compose the columns of the
data set, $D\in\mathbb{R}^{(9\mathrm{x}3)}$. The numbers in the last column of
Table 4 are the squared Mahalanobis distances, $d_{j}$, associated with each
trivariate sample, calculated using (21). Figure 12 provides a quantile-
quantile (QQ) plot for the data set as a graphical method to test for
multivariate normality. The grey shaded region represents a 95% confidence
region with a sample size of $n=9$ for the squared Mahalanobis distance for
trivariate normal data. The squared Mahalanobis distance approximately follows
a chi-square distribution with $q$-degrees of freedom, three degrees of
freedom in this case, $\chi_{q=3}^{2}$ (Johnson and Wichern 2007). In Addition
to the QQ plot, two statistical hypothesis tests were used to check the
assumption of multivariate normally distributed data. An Anderson-Darling and
a Kolmogorov-Smirnov test were used to assess if sample data, $d_{j}$, can be
assumed to approximately originate from a population with a $\chi_{q=3}^{2}$
distribution using a significance value of $\alpha=0.05$. The p-values for the
test were 0.6054 and 0.7185 for the Anderson-Darling and a Kolmogorov-Smirnov
test respectively. Since all the squared Mahalanobis distances from the data
set reside within these confidence bounds, and the p-values were all greater
than the significance value of $\alpha=0.05$, one can assume that the data set
can be adequately represented as trivariate normally distributed.

Figure 12: QQ plot for checking the multivariate Gaussian assumption. The grey
shaded region represents a 95% confidence region with a sample size of $n=9$
for the squared Mahalanobis distance for trivariate normal data. The squared
Mahalanobis distance follows a chi-square distribution with $q$-degrees of
freedom, three degrees of freedom in this case, $\chi_{q=3}^{2}$. The dashed
grey lines represent the lower and upper confidence bounds of the region
representing 2.5% and 97.5% confidence respectively. The red dash-dotted line
represents the theoretical QQ line for $\chi_{q=3}^{2}$. The blue crosses
represent the squared Mahalanobis distance for each sample, $d_{j}$, from the
multivariate data tabulated in Table 4.

Figure 13 is a 3D plot of the uncertainty on a one-sided, CDF based,
trivariate quantile iso-surface with a confidence value of 95% using the data
in Table 4. The first three columns depict the different views which are
normal to the plane corresponding to the XY, XZ, and YZ axis. The fourth
column provides a 3D perspective of the iso-surface. The three rows depict the
various bootstrap confidence interval methods: basic percentile, bias
corrected percentile, and bias corrected and accelerated (BCa) percentile. A
total number of $b=2000$ bootstrap resamples were performed. Significantly
more distortion of the iso-surface incur when utilizing the BCa bootstrap
method compared to the basic percentile method. These distortions appear as
“waviness” as opposed to smoothness in the iso-surfaces. This result is
counter intuitive as the BCa method has precedence of being the more accurate
bootstrap confidence interval method. For the bivariate case the BCa is the
preferred method as it maintains smoothness of the quantile contour while
correcting for bias and skewness. It is not clear whether this distortion is
resolved by increasing the bootstrap sample size, decreasing the mesh
tessellation size, or just limitations of the BCa method for bootstrap
confidence intervals in higher dimensions. Calculating
$\hat{Q}_{0.95}(\tau=0.90)$ requires a significant amount of computation time
in $\mathbb{R}^{3}$ compared to directly calculating the critical points on
$Q(\tau)$. The hat symbol on $\hat{Q}_{\gamma}$ denotes that this is an
estimated subset of $Q(\tau)$, since $Q(\tau)$ is an infinite set of vectors.

Figure 13: 3D plot of the uncertainty on a one-sided, CDF based, trivariate
quantile iso-surface with a confidence value of 95%. The first three columns
depict the different views which are normal to the plane corresponding to the
XY, XZ, and YZ axis. The fourth column provides a 3D perspective of the iso-
surface. The three rows depict the various bootstrap confidence interval
methods: basic percentile, bias corrected percentile, and bias corrected and
accelerated (BCa) percentile.

Figure 14 is a 3D scatter plot of the trivariate data found in Table 4 along
with the critical points using method 1 from section 4.2 and method 2 from
section 4.4 above. Method 1 is the critical point from the iso-surface using
the BCa bootstrap confidence interval on the tessellated CDF. Method 2 is the
direct numerical solution also using the BCa bootstrap confidence interval
method. Both methods used $b=2000$ bootstrap resamples using a nonparametric
bootstrap resampling scheme. The critical point,
$\boldsymbol{v}_{0.95}(0.90)$, from method 2 is strictly greater than the
critical point from method 1, thus a more conservative value. These results
here agree with the numerical simulation study as Method 2 was more
conservative than Method 1 as the desired $p$-value was greater than the
theoretical value of 0.95 for a confidence level of 95%. However, the
difference between the two critical points is significant. The authors have
found that errors start to become significant when calculating the CDF based
quantiles for trivariate distributions with a coarse tessellation mesh size of
0.1 Errors begin to occur as a result of estimating the iso-surface from the
coarse tessellation of the support. Secondly, estimating the critical point
from this coarse iso-surface leads to additional errors. All the errors
accrued in the standardized domain are thus scaled when transforming back the
original domain of the data. For standard deviations greater than unity, this
implies that the errors are multiplied becoming even greater. The difference
in critical point values from the various bootstrap confidence interval
methods were less significant with errors only in the first decimal place.
Therefore, the errors are predominately associated from the differences in the
methodologies. If the critical points are only needed, method 2 provides the
most accurate solution while being over two orders of magnitude faster
(approximately 664x).

Figure 14: 3D scatter plot of the trivariate data found in Table 4 along with
the one-sided confidence bound on the critical points using method 1 from
section 4.2 and method 2 from section 4.4. Method 1 is the critical point from
the iso-surface and Method 2 is the direct numerical solution from Method 2,
both utilizing a BCa bootstrap confidence interval method.

For multi-axis environmental specification in $\mathbb{R}^{3}$, it is
desirable to provide a single conservative point that bounds the majority of
the population with the highest probability density. Thus, at a frequency of
200.24(Hz), a 95% one-sided confidence interval on the critical point of the
multivariate quantile $Q(\tau=0.90)$ with the highest PDF value is
$\boldsymbol{v}_{0.95}(0.90)=[10.0476,\ 18.2621,\ 4.5783](\text{G})$ using
method 2, see Figure 14. Therefore, this critical value bounds the extreme
responses for a proportion of at least 0.90 from a trivariate population with
a given confidence of 95%. The procedure listed above would be repeated for
all subsequent frequencies from the SRS. Coalescing the values from each
confidence interval over all frequencies are then use to construct the final
multi-axis shock environmental specification in $\mathbb{R}^{3}$.

Table 5 provides tabulated results comparing the methodology for confidence
intervals of critical points on multivariate quantiles to univariate tolerance
intervals for each of the three variates represented as three spatial
directions. One can see that the multivariate statistical interval is strictly
greater than each of the univariate tolerance intervals which is necessitated
by the fact that the confidence interval on the multivariate quantile accounts
for each variate but also the correlations between the variates. Utilizing
Algorithm 1 and a one-sided BCa bootstrap confidence interval, one finds that
the uncertainty on multivariate quantile when using concurrent univariate
tolerance intervals is $\tau_{J}=0.76302$ with 95% confidence. As a
consequence of neglecting correlation between the variates, the actual joint
or multivariate quantile probability is less than the desired quantile
probability of $\tau=0.90$. The multivariate methods presented here do not
have these issues as they account for correlations between the variates. The
bounds for the X and Y variates using the Bonferroni interval are less than
the multivariate method. It is not clear that the Bonferroni method would be
strictly less than the multivariate method. The simultaneous tolerance
interval $(\tilde{T}_{\mathrm{simul},\beta=0.90})$ does provide a coverage,
$\beta=\tau$, with 95% confidence but
$F(\tilde{T}_{\mathrm{simul},\beta=0.90})\neq\tau$. The multivariate method on
the other hand provides a coverage, $\beta\geq\tau$ with
$F(\boldsymbol{v}_{0.95}(0.90))=\tau$ with 95% confidence.

Statistical Interval Type | X (G) | Y (G) | Z (G)
---|---|---|---
| Confidence Interval on Critical Point
---
of Multivariate Quantile (Method 2)
10.0476 | 18.2621 | 4.5783
Univariate Tolerance Interval | 9.0081 | 15.9969 | 4.5694
| Bonferroni Simultaneous
---
Tolerance Interval
9.8128 | 17.7368 | 4.8529
Table 5:  Comparison of the proposed multivariate critical point quantile
confidence intervals, univariate tolerance intervals, and the Bonferroni
simultaneous univariate tolerance interval. All are one-sided statistical
intervals using a confidence value of 95% with either a multivariate or
univariate quantile probability value of $\tau=0.90$

##  6 Future Work and Conclusion

Evaluating a multivariate CDF over a tessellated space in $\mathbb{R}^{3}$ can
take several minutes or longer for a single Monte Carlo simulation which
strongly inhibits numerous simulations that are required for constructing
bootstrap confidence intervals. The problem exponentially gets worse with
further increases in dimensionality. One possibility to alleviate the problem
with computing the CDF over a high dimensional space is to empirically
calculate the CDF. Langrené and Warin (2021) have shown that the multivariate
empirical CDF calculation can be fast with a high degree of accuracy. Future
work should also consider the use of adaptive mesh refinement after a coarse
initial estimates of the CDF based multivariate normal quantile has been
calculated. Using an iterative method to refine where the CDF based
multivariate normal quantile are located within $\mathbb{R}^{q}$ will provide
results with higher accuracy, minimizing unnecessary CDF calculations over
areas where the quantile are not located and minimizing memory usage.
Ultimately, CDF based quantile methods consider the tails of the multivariate
distributions while elliptical quantiles from traditional tolerance regions do
not.

Three numerical simulation studies have been provided with all of them
demonstrating great efficacy. The first numerical study was on quantifying the
uncertainty on the multivariate quantile probability, $\tau$, when using
concurrent one-sided univariate quantile probabilities. The second numerical
study was on the constructing confidence internals on CDF based bivariate
normal quantiles, $\hat{Q}_{\gamma}(\tau)$. The third numerical study was on
constructing confidence internals on critical points on CDF based multivariate
normal quantiles, $\hat{\boldsymbol{v}}_{\gamma}(\tau)$, for
$q\in\left\\{2,3,4\right\\}$. The case study provided here demonstrates a
beneficial workflow example for utilizing the presented methodologies and the
current limitations of them. We have shown the inadequacy of univariate
statistical intervals when working with multivariate data. Multi-axis shock
and vibration testing provides a domain where multivariate quantiles are
salient to properly address conservatism when developing environmental test
specifications. The methodologies and algorithms presented here will also help
practitioners using multivariate data with techniques to address the concept
of a “one-sided” statistical interval for tolerance and quantile applications
in a multivariate space.

##  Appendix A Bounds on the Multivariate Quantile Probability

The easiest bound for the multivariate quantile probability is the trivial
solution when all the correlations between the variates are 1, i.e.,
$1=\rho_{{\scriptscriptstyle XY}}=\rho_{{\scriptscriptstyle
XZ}}=\rho_{{\scriptscriptstyle YZ}}$. For this trivial case, since all the
variates are co-linear, the multivariate quantile probability,
$\tau_{{\scriptscriptstyle XYZ}}$, is identical to the individual univariate
quantile probability values,

$\displaystyle\tau_{{\scriptscriptstyle
XYZ}}=\tau_{i}:i\in\left\\{X,Y,Z\right\\}.$ (18)

A lower bound on the multivariate quantile probability for all three variates
can be calculated using the Bonferroni inequality provided by Meeker et al.
(2023), Section D.7.3,

$\displaystyle\tau_{{\scriptscriptstyle XYZ}}\geq
1-\left(\sum_{i=\left\\{X,Y,Z\right\\}}(1-\tau_{i})\right),$ (19)

which occurs when all the variates are maximally, negatively, correlated:
$-1=\rho_{{\scriptscriptstyle XY}}=\rho_{{\scriptscriptstyle
XZ}}=\rho_{{\scriptscriptstyle YZ}}$. In the case where all the variates are
uncorrelated, i.e., $0=\rho_{{\scriptscriptstyle
XY}}=\rho_{{\scriptscriptstyle XZ}}=\rho_{{\scriptscriptstyle YZ}}$, the
multivariate quantile probability can be calculated as,

$\displaystyle\tau_{{\scriptscriptstyle
XYZ}}=\prod_{i=\left\\{X,Y,Z\right\\}}\tau_{i}.$ (20)

##  Appendix B Calculating the Squared Mahalanobis Distances

The proportion of samples that reside inside the tolerance region can be
determined by using the squared Mahalanobis distance. Calculating the squared
Mahalanobis distances, $d_{j}$, for the $j^{\mathrm{th}}$ multivariate sample,
$\boldsymbol{x}_{j}$, can be performed as such:

$\displaystyle
d_{j}=h(\boldsymbol{x}_{j},\bar{\boldsymbol{x}},S):=(\boldsymbol{x}_{j}-\bar{\boldsymbol{x}})S^{-1}(\boldsymbol{x}_{j}-\bar{\boldsymbol{x}})^{T},$
(21)

where $\bar{\boldsymbol{x}}$ denotes the sample mean vector, $S$ denotes the
sample covariance matrix, and the superscript $T$ denotes the transpose
operator. Both $\bar{\boldsymbol{x}}$ and $S$ are plugin estimates for
$\boldsymbol{\mu}$ and $\Sigma$, respectively. One can also calculate the
squared Mahalanobis distances for all the multivariate samples in the data:

$\displaystyle\boldsymbol{d}=h(D,\bar{\boldsymbol{x}},S):=\mathrm{diag}\bigg{(}(D-\bar{\boldsymbol{x}})S^{-1}(D-\bar{\boldsymbol{x}})^{T}\bigg{)},$
(22)

where $\boldsymbol{d}$ is a vector of all the squared Mahalanobis distances
for every multivariate sample. Therefore the $j^{\mathrm{th}}$ element in
$\boldsymbol{d}$ is associated to the squared Mahalanobis distance for the
$j^{\mathrm{th}}$ multivariate sample in $D$.

##  Appendix C Properties of CDF Based Quantile Contours, Iso-Surfaces, and
Hyper Iso-Surfaces

When the multivariate distribution’s parameters are known, one can numerically
estimate the multivariate quantiles by first meshing the support for the
distribution, commonly done through the tessellation of a $q$-dimensional
Euclidean space. The multivariate CDF can then be calculated for any arbitrary
distribution over the tessellated space using most modern computational
languages such as MATLAB or R. For bivariate and trivariate distributions,
adaptive quadrature based on the methods developed by Drezner (1994), Drezner
and Wesolowsky (1990), and Genz (2004) are often employed. However, when there
are four or more dimensions, a quasi-Monte Carlo integration algorithm based
on methods developed by Genz and Bretz (1999) and Genz and Bretz (2002) are
often used. Most multivariate CDF functions also provide an option to control
the absolute error tolerance, typically on the order of 1e-4 or less. After
the CDF values are calculated over the space, one can then extract locations
over the support where approximately equal CDF values are located, thus
creating the CDF based quantile contours for the bivariate case and CDF based
quantile iso-surfaces in the trivariate case.

Figure 15 builds off the example from Figure 3 by showing some CDF based
bivariate normal quantiles, $Q(\tau)$, from
$\tau\in\left\\{0.05,0.1,\ldots,0.90,0.95\right\\}$, shown as a gradient of
colored contours. Figure 16 provides two different plots at three different
quantile probability values, $\tau\in\left\\{0.01,0.90,0.97\right\\}$. The
solid, dashed, and dotted, contours represent correlation coefficients at,
$\rho\in\left\\{0,0.99,-0.99\right\\}$ respectively. At low quantile
probability values, there exists strong differences between the quantile
contours for various correlation coefficients, 16(a). At larger quantile
probability values the differences between the quantiles from various
correlation coefficients become harder to identify, 16(b). When $\tau=0.90$,
it is still easy to identify the case where $\rho=0.99$, as it is very close
to the univariate normal quantiles of $Q(0.90)$. However, in the cases where
$\rho=0$ and $\rho=-0.99$, it becomes much harder to distinguish between the
two contours as there is a large amount of overlap, but they are different.
The grey dashed lines in Figure 16 represent the univariate normal quantiles,
$Q(\tau):\tau\in\left\\{0.01,0.90\right\\}$, which are the asymptotic bounds
for the CDF based multivariate quantiles. Therefore, the CDF based
multivariate quantiles are always strictly greater than or equal to each
variate’s univariate normal quantile.

Figure 15: CDF Based bivariate normal quantiles from
$\tau=\left\\{0.05,0.1,\ldots,0.90,0.95\right\\}$, shown as a gradient of
colored contours. At low quantile probability values, such as $\tau=0.05$,
there exists strong differences between the quantile contours for various
correlation coefficients. At larger quantile probability values, such as
$\tau=0.95$, the differences between the quantiles from various correlation
coefficients become harder to identify.

(a) CDF based quantile contour at a low quantile probability value,
$\tau=0.01$, shown in purple.

(b) CDF based quantile contours at two different large quantile probability
values, $\tau=\left\\{0.90,0.97\right\\}$, shown in blue and gold
respectively.

Figure 16: The solid, dashed, and dotted contours represent correlation
coefficients at, $\rho=\left\\{0,0.99,-0.99\right\\}$ respectively. At low
quantile probability values, there exists strong differences between the
quantile contours for various correlation coefficients. At larger quantile
probability values the differences between the quantiles for correlation
coefficients $\rho=0$ and $\rho=-0.99$ become much harder to identify. At
$\rho=0.99$ the quantile contours are still noticeably different compared to
the quantile contours at $\rho=0$ and $\rho=-0.99$

The figures below all utilize a bivariate normal distribution due to being
easier to visualize and plot. For trivariate normal distributions, the
quantiles are represented as iso-surfaces for constant quantile probability
values. For multivariate distributions in $\mathbb{R}^{q}:q>3$, the quantiles
are hyper iso-surfaces.

##  Appendix D Coverage for Multivariate Distributions

Coverage, in the most general case, for multivariate distribution is taken
from (8) and is defined below,

$\displaystyle\beta=\int\ldots\int_{\Omega}\bigg{(}f(x_{1},\ldots,x_{q})\
dx_{1},\ldots,dx_{q}\bigg{)}.$ (23)

CDF based quantiles in a multivariate space are not intuitive for two main
reasons. The first being ambiguity as there are many (infinitely) number of
domains, $\Omega$, which satisfy (23). Recall the bivariate case from (8), the
domain is $\Omega=\left\\{\boldsymbol{x}=(x_{1},x_{2})\ |\
f(\boldsymbol{x},\bar{\boldsymbol{x}},S)\leq r\right\\}$, which is an
elliptical region. When the domain is $\Omega=\left\\{(x_{1},x_{2})\ |\
\boldsymbol{x}\leq\tilde{T}_{\mathrm{simul},\beta}\right\\}$, the domain is a
rectangular region in $\mathbb{R}^{2}$ which can be seen as the upper right
hand region starting at the intersection of the univariate normal quantiles in
Figure 5, when the distribution parameter’s are known.

The second reason is related to univariate statistics, where the
$\tau^{\mathrm{th}}$ quantile probability also implies coverage, $\beta=\tau$,
since $\beta=\int_{0}^{Q(\tau)}f(x)\,dx$. However, when the domain, $\Omega$,
is changed from the CDF based quantile in (23) to
$\Omega=\left\\{(x_{1},x_{2})\ |\ \boldsymbol{x}\leq Q(\tau)\right\\}$, the
equality now becomes an inequality, $\beta\geq\tau$. Figure 17 provides an
example for a standard normal bivariate distribution with a correlation of
$\rho=0.9$ and a CDF quantile contour, $Q(\tau=0.7)$. Utilizing one hundred
thousand random multivariate samples, the true value can be estimated
accurately, $\hat{\beta}\approxeq\beta$. From Figure 17, it is evidently clear
that $\beta\neq\tau$, but $\beta\geq\tau$. The limiting case is when the
correlation, $\rho$, approaches 1, then $\beta=\tau$,
i.e.,$\lim\limits_{\rho\to 1}\beta=\tau$.

Figure 17: Example of a standard normal bivariate distribution with a
correlation of $\rho=0.9$ and a CDF quantile contour, $Q(\tau=0.7)$, and
$\hat{\beta}=0.7841$.

## References

  * Wald and Wolfowitz (1946) A. Wald and J. Wolfowitz.  Tolerance Limits for a Normal Distribution.  _The Annals of Mathematical Statistics_ , 17(2):208–215, 1946.  ISSN 0003-4851.  URL https://www.jstor.org/stable/2236039.  Publisher: Institute of Mathematical Statistics.
  * Wald (1943) Abraham Wald.  An Extension of Wilks’ Method for Setting Tolerance Limits.  _The Annals of Mathematical Statistics_ , 14(1):45–55, 1943.  ISSN 0003-4851.  URL https://www.jstor.org/stable/2236001.  Publisher: Institute of Mathematical Statistics.
  * Wald (1942) Abraham Wald.  Setting of Tolerance Limits When the Sample is Large.  _The Annals of Mathematical Statistics_ , 13(4):389–399, 1942.  ISSN 0003-4851.  URL https://www.jstor.org/stable/2235841.  Publisher: Institute of Mathematical Statistics.
  * Meeker et al. (2023) William Q. Meeker, Gerald J. Hahn, and Luis A. Escobar.  _Statistical Intervals: A Guide for Practitioners and Researchers, 2nd Edition | Wiley_.  Wiley, March 2023.  ISBN 978-0-471-68717-7.  URL https://www.wiley.com/en-us/Statistical+Intervals%3A+A+Guide+for+Practitioners+and+Researchers%2C+2nd+Edition-p-9780471687177.
  * Chiang and Hsiao (2021) Chieh Chiang and Chin-Fu Hsiao.  Tolerance interval testing for assessing accuracy and precision simultaneously.  _PLoS ONE_ , 16(2):e0246642, February 2021.  ISSN 1932-6203.  10.1371/journal.pone.0246642.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7864420/.
  * Kenett and Zacks (2021) Ron S. Kenett and Shelemyahu Zacks.  _Modern Industrial Statistics: With Applications in R, MINITAB, and JMP_.  John Wiley & Sons, May 2021.  ISBN 978-1-119-71496-5.  URL https://www.wiley.com/en-us/Modern+Industrial+Statistics:+With+Applications+in+R,+MINITAB,+and+JMP,+3rd+Edition-p-9781119714903.  Google-Books-ID: iYEqEAAAQBAJ.
  * Francq et al. (2020) Bernard G. Francq, Marion Berger, and Charles Boachie.  To tolerate or to agree: A tutorial on tolerance intervals in method comparison studies with BivRegBLS R Package.  _Statistics in Medicine_ , 39(28):4334–4349, December 2020.  ISSN 0277-6715.  10.1002/sim.8709.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7756677/.
  * Lucagbo and Mathew (2023) Michael Daniel Lucagbo and Thomas Mathew.  Rectangular tolerance regions and multivariate normal reference regions in laboratory medicine.  _Biometrical Journal_ , 65(3):2100180, 2023.  ISSN 1521-4036.  10.1002/bimj.202100180.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.202100180.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/bimj.202100180.
  * Young and Mathew (2020) Derek S Young and Thomas Mathew.  Nonparametric hyperrectangular tolerance and prediction regions for setting multivariate reference regions in laboratory medicine.  _Statistical Methods in Medical Research_ , 29(12):3569–3585, December 2020.  ISSN 0962-2802.  10.1177/0962280220933910.  URL https://doi.org/10.1177/0962280220933910.  Publisher: SAGE Publications Ltd STM.
  * Chew (1966) Victor Chew.  Confidence, Prediction, and Tolerance Regions for the Multivariate Normal Distribution.  _Journal of the American Statistical Association_ , 61(315):605–617, September 1966.  ISSN 0162-1459.  10.1080/01621459.1966.10480892.  URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1966.10480892.  Publisher: Taylor & Francis _eprint: https://www.tandfonline.com/doi/pdf/10.1080/01621459.1966.10480892.
  * Krishnamoorthy and Mathew (2009) Kalimuthu Krishnamoorthy and Thomas Mathew.  The Multivariate Normal Distribution.  In _Statistical Tolerance Regions: Theory, Applications, and Computation_ , pages 225–247. John Wiley & Sons, Ltd, 2009.  ISBN 978-0-470-47390-0.  10.1002/9780470473900.ch9.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9780470473900.ch9.  Section: 9 _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470473900.ch9.
  * Bornkamp (2018) Björn Bornkamp.  Calculating quantiles of noisy distribution functions using local linear regressions.  _Computational Statistics_ , 33(1):487–501, March 2018.  ISSN 1613-9658.  10.1007/s00180-017-0736-0.  URL https://doi.org/10.1007/s00180-017-0736-0.
  * Chen and Welsh (2002) L. A. Chen and A. H. Welsh.  Distribution-Function-Based Bivariate Quantiles.  _Journal of Multivariate Analysis_ , 83(1):208–231, October 2002.  ISSN 0047-259X.  10.1006/jmva.2001.2043.  URL https://www.sciencedirect.com/science/article/pii/S0047259X01920433.
  * Klein and Kneib (2020) Nadja Klein and Thomas Kneib.  Directional bivariate quantiles: a robust approach based on the cumulative distribution function.  _AStA Advances in Statistical Analysis_ , 104(2):225–260, June 2020.  ISSN 1863-818X.  10.1007/s10182-019-00355-3.  URL https://doi.org/10.1007/s10182-019-00355-3.
  * Efron (1979) B. Efron.  Bootstrap Methods: Another Look at the Jackknife.  _The Annals of Statistics_ , 7(1):1–26, 1979\.  ISSN 0090-5364.  URL https://www.jstor.org/stable/2958830.  Publisher: Institute of Mathematical Statistics.
  * Efron and Tibshirani (1986) B. Efron and R. Tibshirani.  Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy.  _Statistical Science_ , 1(1):54–75, 1986.  ISSN 0883-4237.  URL https://www.jstor.org/stable/2245500.  Publisher: Institute of Mathematical Statistics.
  * Efron and Hastie (2016) Bradley Efron and Trevor Hastie.  The Jackknife and the Bootstrap.  In _Computer Age Statistical Inference: Algorithms, Evidence, and Data Science_ , Institute of Mathematical Statistics Monographs, pages 155–180. Cambridge University Press, Cambridge, 2016.  ISBN 978-1-107-14989-2.  10.1017/CBO9781316576533.011.  URL https://www.cambridge.org/core/books/computer-age-statistical-inference/jackknife-and-the-bootstrap/900BA9154E0F7F487FDE155336A99BF1.
  * Langrené and Warin (2021) Nicolas Langrené and Xavier Warin.  Fast multivariate empirical cumulative distribution function with connection to kernel density estimation.  _Computational Statistics & Data Analysis_, 162:107267, October 2021.  ISSN 0167-9473.  10.1016/j.csda.2021.107267.  URL https://www.sciencedirect.com/science/article/pii/S0167947321001018.
  * Davison and Hinkley (1997) A. C. Davison and D. V. Hinkley.  Bootstrap Confidence Intervals.  In _Bootstrap Methods and their Application_ , Cambridge Series in Statistical and Probabilistic Mathematics, pages 191–255. Cambridge University Press, Cambridge, 1997.  ISBN 978-0-521-57471-6.  10.1017/CBO9780511802843.006.  URL https://www.cambridge.org/core/books/bootstrap-methods-and-their-application/confidence-intervals/499FDF2036D956E7B7E7FB7F25D08B13.
  * Carpenter and Bithell (2000) James Carpenter and John Bithell.  Bootstrap confidence intervals: when, which, what? A practical guide for medical statisticians.  _Statistics in Medicine_ , 19(9):1141–1164, 2000\.  ISSN 1097-0258.  10.1002/(SICI)1097-0258(20000515)19:9¡1141::AID-SIM479¿3.0.CO;2-F.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0258%2820000515%2919%3A9%3C1141%3A%3AAID-SIM479%3E3.0.CO%3B2-F.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/%28SICI%291097-0258%2820000515%2919%3A9%3C1141%3A%3AAID-SIM479%3E3.0.CO%3B2-F.
  * Lewandowski et al. (2009) Daniel Lewandowski, Dorota Kurowicka, and Harry Joe.  Generating random correlation matrices based on vines and extended onion method.  _Journal of Multivariate Analysis_ , 100(9):1989–2001, October 2009.  ISSN 0047-259X.  10.1016/j.jmva.2009.04.008.  URL https://www.sciencedirect.com/science/article/pii/S0047259X09000876.
  * noa (a) TMIL-STD-810 TEST METHOD STANDARD: ENVIRONMENTAL ENGINEERING CONSIDERATIONS AND LABORATORY TESTS, a.
  * noa (b) JESD22S Solid State Device Packaging Standards B103 (Vibration, Variable Frequency), b.  URL https://www.desolutions.com/testing-services/test-standards/jesd22/.
  * noa (c) NAVMAT P-9492 NAVY MANUFACTURING SCREENING PROGRAM DECREASE CORPORATE COSTS INCREASE FLEET READINESS, c.
  * Whiteman and Berman (2002) Wayne E. Whiteman and Morris S. Berman.  Fatigue failure results for multi-axial versus uniaxial stress screen vibration testing.  _Shock and Vibration_ , 9(6):319–328, January 2002.  ISSN 1070-9622.  URL https://content.iospress.com/articles/shock-and-vibration/sav00185.  Publisher: IOS Press.
  * Sisemore et al. (2020) Carl Sisemore, Vit Babuska, and Robert-X. Flores.  Multi-Axis Resonant Plate Shock Testing Evaluation and Test Specification Development.  Technical Report SAND-2020-10224, Sandia National Lab. (SNL-NM), Albuquerque, NM (United States), September 2020.  URL https://www.osti.gov/biblio/1670873.
  * Nelson (2017) Garrett Dean Nelson.  A Systematic Evaluation of Test Specification Derivation Methods for Multi-Axis Vibration Testing.  Technical Report SAND2017-11538C, Sandia National Lab. (SNL-NM), Albuquerque, NM (United States), October 2017.  URL https://www.osti.gov/biblio/1480216.  ISSN: 2191-5644.
  * Lalanne (2014) Christian Lalanne.  _Mechanical Vibration and Shock Analysis, Volume 5, Specification Development, 3rd Edition | Wiley_.  Wiley, April 2014.  ISBN 978-1-118-93123-3.  URL https://www.wiley.com/en-ie/Mechanical+Vibration+and+Shock+Analysis%2C+Volume+5%2C+Specification+Development%2C+3rd+Edition-p-9781118931233.
  * Jacobs-Omalley (2014) Laura Diane Jacobs-Omalley.  Methodology for Defining Multi-Axis Vibration Specifications.  Technical Report SAND2014-18809C, Sandia National Lab. (SNL-NM), Albuquerque, NM (United States), October 2014.  URL https://www.osti.gov/biblio/1315349.
  * Nelson and Jacobs-Omalley (2015) Garrett Dean Nelson and Laura Diane Jacobs-Omalley.  Comparison of Multi-Axis and Single Axis Testing on Plate Structures.  Technical Report SAND2015-1363C, Sandia National Lab. (SNL-NM), Albuquerque, NM (United States), February 2015.  URL https://www.osti.gov/biblio/1240323.
  * Johnson and Wichern (2007) Richard A. Johnson and Dean W. Wichern.  _Applied Multivariate Statistical Analysis_.  Pearson, Upper Saddle River, N.J, 6th edition edition, April 2007.  ISBN 978-0-13-187715-3.
  * Drezner (1994) Zvi Drezner.  Computation of the Trivariate Normal Integral.  _Mathematics of Computation_ , 62(205):289–294, 1994.  ISSN 0025-5718.  10.2307/2153409.  URL https://www.jstor.org/stable/2153409.  Publisher: American Mathematical Society.
  * Drezner and Wesolowsky (1990) Zvi Drezner and G. O. Wesolowsky.  On the computation of the bivariate normal integral.  _Journal of Statistical Computation and Simulation_ , 35(1-2):101–107, March 1990.  ISSN 0094-9655.  10.1080/00949659008811236.  URL https://doi.org/10.1080/00949659008811236.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/00949659008811236.
  * Genz (2004) Alan Genz.  Numerical computation of rectangular bivariate and trivariate normal and t probabilities.  _Statistics and Computing_ , 14(3):251–260, August 2004.  ISSN 1573-1375.  10.1023/B:STCO.0000035304.20635.31.  URL https://doi.org/10.1023/B:STCO.0000035304.20635.31.
  * Genz and Bretz (1999) Alan Genz and Frank Bretz.  Numerical computation of multivariate t-probabilities with application to power calculation of multiple contrasts.  _Journal of Statistical Computation and Simulation_ , 63(4):103–117, July 1999.  ISSN 0094-9655.  10.1080/00949659908811962.  URL https://doi.org/10.1080/00949659908811962.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/00949659908811962.
  * Genz and Bretz (2002) Alan Genz and Frank Bretz.  Comparison of Methods for the Computation of Multivariate t Probabilities.  _Journal of Computational and Graphical Statistics_ , 11(4):950–971, December 2002.  ISSN 1061-8600.  10.1198/106186002394.  URL https://doi.org/10.1198/106186002394.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1198/106186002394.