\section{Introduction}
Semiparametric copulas are useful tools in multivariate data analysis. They are used for modelling a multivariate distribution whose dependence structure is induced by a known copula and whose marginal distributions are estimated; see, for example, Sancetta and Satchell (2004), Segers et al. (2014) and the references therein.  We point out to the interesting  work of Rosen and Thompson (2015) who proposed a semiparametric methodology for modeling a multivariate distribution whose dependence structure is induced by a Gaussian copula and whose marginal distributions are estimated nonparametrically via mixtures of B-spline densities.  The authors take a Bayesian approach, using Markov chain Monte Carlo methods for inference.

In the present paper, a Bayesian Semiparametric copula approach based on the Dirichlet process and the Gaussian copula is proposed to model the underlying multivariate distribution $F_{true}$. In addition, recognizing that  many recent applications of research are developed based on the assumption of multivariate normality (Fernandez, 2010 and Zhu et al., 2014), a test  to assess this assumption is developed. Recent procedures tackling this problem can be found in Kim and Park (2018), Madukaife and Okafor (2018), Henze and Visagie (2019) and Al-Labadi et al. (2019a). We highlight that, while most available works in the area of the hypothesis testing using copula approaches are related to assess  independence  (Genest and R\'{e}millard, 2004; Kojadinovic and Holmes, 2009; Medovikov, 2016; Belalia et al., 2017), the proposed test is Bayesian and considers modeling the dependence structure and the marginal behaviors of the data separately to assess the multivariate normality assumption. Briefly, all univariate marginal distributions of $F_{true}$ are assumed to have the Dirichlet process to define posterior-based and prior-based models of $F_{true}$.  A Gaussian copula model is then utilized to induce the dependence structure. The test follows by comparing the concentration of the posterior-based model to the concentration of the prior-based model about the family of multivariate normal distributions (hypothesized model) via the so-called relative belief ratio.  In this comparison, a Bayesian counterpart of the Energy distance is developed. The proposed test is easy to implement with a powerful performance and allows to state evidence for or against the null hypothesis. Also, unlike the test presented in Al-Labadi et al. (2019a),  which is restricted to assess the family of multivariate normal distributions for the hypothesized model, the proposed test can be extended to assess all families of multivariate distributions (model checking problems).

The rest of the paper is structured as follows. A relevant background containing some definitions and generic properties are reviewed in Section 2. In Section 3, a Bayesian semiparametric Gaussian copula approach based on the Dirichlet process is proposed for modeling multivariate distributions. The choice of the hyperparameter of the Dirichlet process and the estimation method of the parameter of the Gaussian copula are discussed in Section. In Section 5, a Bayesian multivariate normality (MVN) test based on the proposed approach and the Energy distance is developed. The main steps of a computational algorithm to implement the MVN test are outlined in Section 6. The performance of the approach and its application to the MVN test is clarified through some simulation studies and a real data example in Section 7. The results show that the proposed test works well in all covered cases and it is very powerful. Finally, Section 8 concludes the paper with a summary of the results. Some notations related to the Section 7 are given in the Appendix.

\section{Relevant background}
\subsection{Copula-based Model}

In multivariate analysis, copula models are introduced by Sklar (1959) as a common tool to model multivariate distributions. Following Nelsen (2006),
an $m$-dimensional copula ($m$-copula) is a nondecreasing and right continuous function $C$ from $[0,1]^{m}$ into $[0,1]$ such that, for every
$\mathbf{u}=(u_{1},\ldots,u_{m})\in [0,1]^{m}$
\begin{enumerate}[nolistsep,label=(\roman*),leftmargin=\parindent,align=left,labelwidth=\parindent,labelsep=0pt]
\item
$C(u_{i},\ldots,u_{i-1},0,u_{i+1},\ldots,u_{m})=0$, for $i=1,\ldots,m$ ($C$ is grounded).
\item
$C(1,\ldots,1,u_{i},1,\ldots,1)=u_{i}$, for $i=1,\ldots,m$ ($C$ has margins).
\item
\hspace{0.1mm} For every $\mathbf{v}=(v_{1},\ldots,v_{m})\in [0,1]^{m}$ such that $u_{i}\leq v_{i}$ for all $i$, the $C$-volume
$V\left([\mathbf{u},\mathbf{v}]\right)\geq 0$ ($C$ is $m$-increasing), where $m$-cube $[0,1]^{m}$ is $m$ product of $[0,1]$ and $[\mathbf{u},\mathbf{v}]$ is the $m$-box $[u_{1},v_{1}]\times\cdots\times[u_{m},v_{m}]$.
\end{enumerate}

From (iii), it is obvious that every $m$-copula $C$ is nondecreasing in each variable and satisfies in the Lipschitz condition. That is, for every point $\mathbf{u}, \mathbf{v}\in [0,1]^{m}$,
\begin{equation}\label{Lipschitz}
\big| C(\mathbf{u})-C(\mathbf{v}) \big| \leq \sum_{i=1}^{m}|u_{i}-v_{i}|.
\end{equation}

\noindent Hence, any $m$-copula $C$ is uniformly continuous on $[0,1]^{m}$.

The following key theorem of Sklar (1959) illustrates the role of the $m$-copulas to model the multivariate distribution functions through their univariate margins.

\begin{theorem}\label{sklar's Thm}
\textbf{(Sklar's theorem)} Let $F$ be an $m$-variate distribution function with marginal distribution functions $F_{1},\ldots,F_{m}$. Then there exists an $m$-copula $C$ such that for all $\mathbf{t}\in \overline{\mathbb{R}}^{m}$
\begin{equation}\label{sklar}
F(t_{1},\ldots, t_{m})=C\left(F_{1}(t_{1}),\ldots, F_{m}(t_{m})\right),
\end{equation}
where $\overline{\mathbb{R}}^{m}$ is $m$ product of the extended real line $[-\infty,\infty]$. If $F_{1},\ldots,F_{m}$ are all continuous, then $C$ is unique and can be written as
 \begin{equation*}
 C(\mathbf{u})=F\left(F^{-1}_{1}(u_{1}),\ldots, F_{m}^{-1}(u_{m})\right),
 \end{equation*}
for any $\mathbf{u}\in[0,1]^{m}$ where $F^{-1}_{i}(u_{i})=\inf\lbrace t\in \mathbb{R}|\, F_{i}(t)\geq u_{i}\rbrace$, otherwise; $C$ is uniquely determined on
$Ran (F_{1})\times \cdots \times Ran (F_{m})$, where $Ran (F_{i})$ denotes the range of the distribution function $F_{i}$ for $i=1,\ldots,m$. Conversely, if $C$ is an $m$-copula and $F_{1},\ldots,F_{m}$ are distribution functions, then $F$, defined by \eqref{sklar}, is an $m$-variate distribution function with marginal distribution functions $F_{1},\ldots,F_{m}$.
\end{theorem}

In general, the dependence structure of multivariate distributions is modeled by the copula. For this purpose, some families of copulas have been developed. See, for example, the work of Joe (1997), Trivedi and Zimmer (2005), and Nelsen (2006). One instance of special interest is the family of \emph{Gaussian copulas}. Beside that it satisfies both Fr\'{e}chet-Hoeffding lower and upper bounds (Nelsen 2006, Theorem 2.10.2), it has only one dependence parameter restricted to the symmetric interval $[-1,1]$, which makes it simple to apply. Formally, let $\Phi^{-1}$ be the inverse of the cumulative distribution function (cdf) of the univariate standard normal distribution and $\Phi_{R}$ be the cdf of the multivariate normal distribution with zero mean vector $\mathbf{0}_{m}$ and correlation matrix $R=(r_{ij})$ for $1\leq i,j\leq m$, then the family of Gaussian copulas is defined by
\begin{equation*}\label{gaussian-cop}
\left\lbrace C_{R}(\mathbf{u}):=\Phi_{R}( \Phi^{-1}(u_{1}),\ldots, \Phi^{-1}(u_{m}))|\hspace{.2cm} \, r_{ij}\in[-1,1]\right\rbrace,
\end{equation*}
\noindent where $\mathbf{u}\in [0,1]^{m}$. Following Chen et al. (2006), assume that, for any $\mathbf{u}\in [0,1]^{m}$, $C(\mathbf{u})=C_{R}(\mathbf{u})$, then the multivariate distribution $F$ is of a \emph{semiparametric Gaussian copula model} $F(t_{1},\ldots,t_{m})=C_{R}(F_{1}(t_{1}),\ldots, F_{m}(t_{m}))$, with unknown parameter $R$ and unknown marginal cdf $F_{i}$, for $i=1\ldots m$. A detailed discussion about the semiparametric Gaussian copula model will be presented in Section 3 based on using the Dirichlet process.

The following algorithm shows the steps of generating a sample of random vectors from an $m$-variate distribution $F$ with marginal cdf's $F_{1},\ldots,F_{m}$ using a Gaussian copula model with correlation matrix $R$ .
\begin{algorithm}
\caption{Generating a sample from an $m$-variate distribution $F$ using Gaussian copula model}\label{alg1
\begin{enumerate}
\item Generate a random vector $(y_{1},\ldots,y_{m})$ from an $m$-variate normal distribution with mean vector $\mathbf{0}_{m}$ and correlation matrix of $R$. \textbf{\smallskip}
\item Compute $u_{i}=\Phi(y_{i})$, for $i=1,\ldots,m$. \textbf{\smallskip}
\item Compute vector $(x_{1},\ldots,x_{m})$ such that $x_{i}=F^{-1}_{i}(u_{i})$, for $i=1,\ldots,m$.
\item Repeat steps 1-3 for $k$ times to generate a sample of size $k$ from distribution $F$.\textbf{\smallskip}
\end{enumerate}
\end{algorithm}
\subsection{Dirichlet Process}
The Dirichlet process prior, introduced by Ferguson (1973), is the
most commonly used prior in Bayesian nonparametric inferences. A remarkable collection of nonparametric inferences have been devoted to this prior. Here we only present the most relevant
definitions and properties of this prior. Consider a space
$\mathfrak{X}$ with a $\sigma$-algebra $\mathcal{A}$ of subsets of
$\mathfrak{X}$, let $H$ be a fixed probability measure on $(\mathfrak{X}%
,\mathcal{A}),$ called the \emph{base measure}, and $a$ be a positive number,
called the \emph{concentration parameter}. A random probability measure
$P=\left\{  P(A):A\in\mathcal{A}\right\}  $ is called a Dirichlet process on
$(\mathfrak{X},\mathcal{A})$ with parameters $a$ and $H,$ denoted by $P\sim
{DP}(a,H),$ if for every measurable partition $A_{1},\ldots,A_{k}$ of
$\mathfrak{X} $ with $k\geq2\mathfrak{,}$ the joint distribution of the vector
$\left(  P(A_{1}),\ldots\,P(A_{k})\right)$ has the Dirichlet distribution with parameter
$aH(A_{1}),\ldots,$ $aH(A_{k})$. Also, it is assumed that
$H(A_{j})=0$ implies $P(A_{j})=0$ with probability one. Consequently, for any
$A\in\mathcal{A}$, $P(A)\sim$ beta$(aH(A),a(1-H(A)))$,
${E}(P(A))=H(A)\ $and ${Var}(P(A))=H(A)(1-H(A))/(1+a).$ Accordingly, the base measure $H$
plays the role of the center of $P$ while the concentration parameter $a$ controls the variation of
 $P$ around  $H$. One of the most well-known
properties of the Dirichlet process is the conjugacy property. That is, when the sample $x=(x_{1},\ldots,x_{n})$
is drawn from $P\sim DP(a,H)$, the posterior distribution of $P$ given $x$,
denoted by $P^{\ast}$, is also a Dirichlet process with
concentration parameter $a+n$ and base measure
\begin{equation}\label{pos base measure}
H^{\ast}=a(a+n)^{-1}H+n(a+n)^{-1}F_{n},
\end{equation}
where $F_{n}$ denotes the empirical cumulative distribution function (cdf) of the sample
$x$. Note that, $H^{\ast}$ is a convex combination
of the base measure $H$ and the empirical cdf $F_{n}$. Therefore, $H^{\ast}\rightarrow H$ as
$a\rightarrow\infty$ while $H^{\ast}\rightarrow F_{n}$ as $a\rightarrow0.$ A guideline
about choosing the hyperparameters $a$ and $H$ will be covered in Section 4. Following Ferguson (1973), $P\sim{DP}(a,H)\ $can be represented as
\begin{equation}
P=\sum_{i=1}^{\infty}L^{-1}(\Gamma_{i}){\delta_{Y_{i}}/}\sum_{i=1}^{\infty
}{{L^{-1}(\Gamma_{i})}}, \label{series-dp}%
\end{equation}
where $\Gamma_{i}=E_{1}+\cdots+E_{i}$ with $E_{i}\overset{i.i.d.}{\sim}%
$\ exponential$(1),Y_{i}\overset{i.i.d.}{\sim}H$ independent of the
$\Gamma_{i},L^{-1}(y)=\inf\{x>0:L(x)\geq y\}$ with $L(x)=a\int_{x}^{\infty
}t^{-1}e^{-t}dt,x>0,$ and ${\delta_{a}}$ the Dirac delta measure. The series representation
(\ref{series-dp}) implies that the Dirichlet process is a discrete probability
measure even for the cases with an absolutely continuous base measure $H$. Note
that, by imposing the weak topology, the support of the Dirichlet process could
be quite large. Recognizing the complexity when working with \eqref{series-dp}, Zarepour and Al-Labadi (2012) proposed the following finite representation as an efficient method to simulate the Dirichlet process. They showed that the Dirichlet process $P\sim DP(a,H)$ can be approximated by
\begin{equation}\label{approx of DP}
P_{N}=\sum_{i=1}^{N}J_{i,N}\delta_{Y_{i}},
\end{equation}
with the
monotonically decreasing weights
$J_{i,N}=\frac{G_{a/N}^{-1}(\frac{\Gamma_{i}}{\Gamma_{N+1}})}{\sum_{j=1}^{N}G_{a/N}^{-1}(\frac{\Gamma_{i}}{\Gamma_{N+1}})},$
where $\Gamma_{i}$ and $Y_{i}$ are defined as before, $N$ is a positive large
integer and $G_{a/N}$ denotes the complement-cdf of the $\text{gamma}(a/N,1)$ distribution.
Note that, $G^{-1}_{a/N}(p)$ is the $(1-p)$-th quantile of the $\text{gamma}(a/N,1)$ distribution.
The following
algorithm describes how the approximation \eqref{approx of DP} can be used to generate
a sample from $DP(a,H)$.

\begin{algorithm}
\caption{Approximately generating a sample from $DP(a,H)$}
\begin{enumerate}
\item Fix a large positive integer $N$ and generate i.i.d. $Y_{i}\sim H$ for $i=1,\ldots,N$.

\item For $i=1,\ldots,N+1$, generate i.i.d. $E_{i}$ from the exponential distribution with rate 1, independent
of $\left(Y_{i}\right)_{1\leq i\leq N}$ and put $\Gamma_{i}=E_{1}+\cdots+E_{i}$.
\item Compute $G_{a/N}^{-1}\left(  {\Gamma_{i}}/{\Gamma_{N+1}}\right)  $ for
$i=1,\ldots,N$ and return $P_{N}.$
\end{enumerate}
\end{algorithm}
The Dirichlet process can also be obtained from the following finite mixture models developed by Ishwaran and Zarepour (2002). Let $P_{N}$ has the from given \eqref{approx of DP} with $(J_{i,N})_{1\leq i\leq N}\sim$ Dirichlet$(a/N,\ldots,a/N)$. Then $E_{P_{N}}(g)\rightarrow E_{P}(g)$ in distribution as $N\rightarrow\infty$, for any measurable function $g: \mathbb{R}\rightarrow \mathbb{R}$ with $\int_{\mathbb{R}}|g(x)|\, H(dx)<\infty$ and $P\sim DP(a,H)$. In particular, $(P_{N})_{N\geq 1}$ converges in distribution to $P$, where $P_{N}$ and $P$ are random values in the space $M_{1}(\mathbb{R})$ of probability measures on $\mathbb{R}$ endowed with the topology of weak convergence. To generate $(J_{i,N})_{1\leq i\leq N}$ put $J_{i,N}=G_{i,N}/\sum_{i=1}^{N}G_{i,N}$, where $(G_{i,N})_{1\leq i\leq N}$ is a sequence of i.i.d. gamma$(a/N, 1)$ random variables independent of $(Y_{i})_{1\leq i\leq N}$. This form of approximation leads to some results in Section 5.
\subsection{Relative Belief Inferences}

The relative belief ratio, developed by Evans (2015), becomes a widespread measure of statistical evidence. See, for example, the work of Al-Labadi and Evans (2018), Al-Labadi et al. (2017, 2018), and Al-Labadi et al. (2019a, 2019b) for implementation of the relative belief ratio on different stimulating univariate hypothesis testing problems. In details, let $\{f_{\theta}:\theta\in\Theta\}$ be a collection of densities on a sample space $\mathfrak{X}$ and let $\pi$ be a prior on the parameter space $\Theta$. Note that
the densities may represent discrete or continuous probability measures but they are
all with respect to the same support measure  $d\theta$. After
observing the data $x,$ the posterior distribution of $\theta$, denoted by $\pi(\theta\,|\,x)$, is a revised prior and is given by the
density $\pi(\theta\,|\,x)=\pi(\theta)f_{\theta}(x)/m(x)$, where $m(x)=\int
_{\Theta}\pi(\theta)f_{\theta}(x)\,d\theta$ is the prior predictive density of
$x.$  For a parameter of interest $\psi=\Psi(\theta),$ let $\Pi_{\Psi}$ be
the marginal prior probability measure and $\Pi_{\Psi}(\cdot|\,x)$ be
the marginal posterior probability measure. It is assumed that
$\Psi$ satisfies regularity conditions
so that the prior density $\pi_{\Psi}$ and the posterior density
$\pi_{\Psi}(\cdot\,|\,x)$ of $\psi$ exist with respect to some support measure on the range space for $\Psi$
. The relative belief ratio for a value
$\psi$ is then defined by $RB_{\Psi}(\psi\,|\,x)=\lim_{\delta\rightarrow0}%
\Pi_{\Psi}(N_{\delta}(\psi\,)|\,x)/\Pi_{\Psi}(N_{\delta}(\psi\,)),$ where
$N_{\delta}(\psi\,)$ is a sequence of neighborhoods of $\psi$ converging
nicely to $\psi$ as $\delta\rightarrow0$ (Evans, 2015). When $\pi_{\Psi}$ and  $\pi_{\Psi}(\cdot\,|\,x)$ are continuous at $\psi,$ the relative belief ratio is defined by
\begin{equation*}
RB_{\Psi}(\psi\,|\,x)=\pi_{\Psi}(\psi\,|\,x)/\pi_{\Psi}(\psi), \label{relbel}%
\end{equation*}
the ratio of the posterior density to the prior density at $\psi.$  Therefore,
$RB_{\Psi}(\psi\,|\,x)$ measures the change in the belief of $\psi$ being the true value from a \textit{priori} to a \textit{posteriori}.

Since $RB_{\Psi}(\psi\,|\,x)$ is a measure of the evidence that $\psi$ is the true value, if $RB_{\Psi}(\psi\,|\,x)$ $>1$, then the probability of the $\psi$ being the true value from a priori to a posteriori is increased, consequently there is evidence based on the data that $\psi$ is the true value. If $RB_{\Psi}(\psi\,|\,x)<1$, then the probability of the $\psi$ being the true value from a priori to a posteriori is decreased. Accordingly, there is evidence against based on the data that $\psi$ being the true value. For the case $RB_{\Psi}(\psi\,|\,x)=1$ there is no
evidence either way.

Obviously, $RB_{\Psi}(\psi_{0}\,|\,x)$ measures the evidence of the hypothesis $\mathcal{H}_{0}:\Psi(\theta)=\psi_{0}$. Large values of $RB_{\Psi}(\psi_{0}\,|\,x)=c$ provides
 strong evidence in favor of $\psi_{0}$. However, there may also exist other
values of $\psi$ that had even larger increases. Thus, it is also necessary, however, to calibrate whether this is strong or weak evidence for
or against $\mathcal{H}_{0}.$ A typical
calibration of $RB_{\Psi}(\psi_{0}\,|\,x)$  is given by the  \textit{strength}
\begin{equation}
\Pi_{\Psi}\left[RB_{\Psi}(\psi\,|\,x)\leq RB_{\Psi}(\psi_{0}\,|\,x)\,|\,x\right].
\label{strength}%
\end{equation}
The value in \eqref{strength} indicates that the posterior probability that the true value of $\psi$ has a relative
belief ratio no greater than that of the hypothesized value $\psi_{0}.$ Noticeably, (\ref{strength}) is not a p-value as it has a very different
interpretation. When $RB_{\Psi}(\psi_{0}\,|\,x)<1$, there is evidence
against $\psi_{0},$ then a small value of (\ref{strength}) indicates
 strong evidence against $\psi_{0}$. On the other hand, a large value for \eqref{strength}    indicates   weak evidence against $\psi_{0}$.
Similarly, when $RB_{\Psi}(\psi_{0}\,|\,x)>1$, there is  evidence in favor
of $\psi_{0},$ then a small value of (\ref{strength}) indicates  weak
evidence in favor of $\psi_{0}$, while a large value of \eqref{strength} indicates
 strong evidence in favor of $\psi_{0}$.
\subsection{Energy Distance}
The \emph{Energy distance}, presented by Sz\'{e}kely (2003), is an appropriate tool to determine the equality of distributions. In general, the Energy distance between two $m$-variate distribution function $F$ and $G$ is defined by
\begin{equation}\label{E-distance-1}
d_{\mathcal{E}}(F,G)=2E\|\mathbf{X}-\mathbf{Y}\|-E\|\mathbf{X}-\mathbf{X}^{\prime}\|-E\|\mathbf{Y}-\mathbf{Y}^{\prime}\|,
\end{equation}
where $\mathbf{X},\mathbf{X}^{\prime}\overset{i.i.d.}{\sim} F$, $\mathbf{Y},\mathbf{Y}^{\prime}\overset{i.i.d.}{\sim} G$ and $\|\mathbf{a}\|=\sqrt{\mathbf{a}^{T}\mathbf{a}}$ denotes  Euclidean norm of vector $\mathbf{a}=(a_{1},\ldots,a_{m})$. Sz\'{e}kely and Rizzo (2013) showed that $d_{\mathcal{E}}(F,G)\geq 0$ such that equality holds if and only if $F=G$. Note that, from (Sz\'{e}kely, 2003), the Energy distance \eqref{E-distance-1} is rotation invariant. This property makes it appropriate for testing goodness-of-fit problems in higher dimensions. Specifically, let  $G$ be the hypothesized distribution and $\mathbf{x}_{m\times n}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$ be the observed sample from $F$. Then, the one sample Energy distance corresponding to \eqref{E-distance-1} is defined by
\begin{equation}\label{E-distance-2}
d_{\mathcal{E},n}(F,G)=\frac{2}{n}\sum_{i=1}^{n}E||\mathbf{x}_{i}-\mathbf{Y}||-\frac{1}{n^{2}}\sum_{\ell,m=1}^{n}||\mathbf{x}_{\ell}-\mathbf{x}_{m}||-E||\mathbf{Y}-\mathbf{Y}^{\prime}||,
\end{equation}
where $\mathbf{x}_{i}\in\mathbb{R}^{m}$, for $i=1,\ldots,n$, and the expectations are taken with respect to the distribution $G$. The special important case occurs when $G$ is a multivariate normal distribution where the $\mathsf{R}$ package \textbf{energy} is usually used for implementing \eqref{E-distance-2}. For further discussion about Energy distance consult Sz\'{e}kely (2003) and Sz\'{e}kely and Rizzo (2013).

\section{A Bayesian Semiparametric Gaussian Copula Approach for Modeling Multivariate Distributions}
In this section, we propose a Bayesian semiparametric copula approach based on the Gaussian copula as a flexible model for modeling multivariate distributions. For a brevity, we refer to this procedure as the BSPGC (Bayesian semiparametric Gaussian copula) approach. Specifically, let $\mathbf{x}_{m\times n}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$ be a sample of size $n$ from an unknown $m$-variate distribution $F_{true}$ with maginal cdf's $F_{1},\ldots,F_{m}$. Note that, the subscript $m\times n$ may be omitted whenever it is clear in the context.
To model $F_{true}$ based on the BSPGC approach we use the prior $F_{i}\sim DP(a,H_{i})$, where $H_{i}$ is the $i$-th marginal cdf of a given $m$-variate distribution $H$. So, by \eqref{pos base measure}, for a given choice of $a$, $F^{\ast}_{i}=F_{i}|\mathbf{x}_{i}\sim DP(a+n, H_{i}^{\ast})$ for $i=1,\ldots,m$.  Consider the joint cdf $H^{\ast}$ corresponding to marginal cdf's $H_{1}^{\ast},\ldots,H_{m}^{\ast}$ with correlation matrix $R^{\ast}$. Then, the \emph{posterior-based model} is defined by
\begin{equation}\label{pos-BSPGC}
F^{\ast}(t_{1},\ldots,t_{m})=C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m})).
\end{equation}
The next lemma shows that $F^{\ast}$ approaches to the true distribution $F_{true}$ when the sample size increases.
\begin{lemma}\label{post as n to inf}
Let $\mathbf{x}_{m\times n}$ be a sample from $m$-variate distribution function $F_{true}$ with unknown marginal cdf's $F_{1},\ldots,F_{m}$. Assume that $F^{\ast}_{i}\sim DP(a+n, H_{i}^{\ast})$, for $ i=1,\ldots,m$. For any $\mathbf{t}=(t_{1},\ldots,t_{m})\in \mathbb{R}^{m}$, $C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m}))\xrightarrow{a.s.} F_{true}(\mathbf{t})$ as $n\rightarrow\infty$.
\end{lemma}
\begin{proof}
For any $\mathbf{t}\in\mathbb{R}^{m}$, the triangle inequality implies
\begingroup\makeatletter\def\f@size{9}\check@mathfonts
  \begin{align}\label{triangular}
\big|C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m}))-F_{true}(\mathbf{t})\big|&\leq\big|C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m}))-H^{\ast}(\mathbf{t})\big|\nonumber\\
&+ \big|H^{\ast}(\mathbf{t})-F_{true}(\mathbf{t})\big|=I_{1}+I_{2}.
\end{align}
\endgroup
From \eqref{pos base measure}, for any $\mathbf{t}\in\mathbb{R}^{m}$, $H^{\ast}(\mathbf{t})\xrightarrow{a.s.} F_{true}(\mathbf{t})$ as $n\rightarrow\infty$. Hence, the continuous mapping theorem implies $I_{2}\xrightarrow{a.s.} 0$ as $n\rightarrow\infty$. On the other hand, from Sklar's theorem \ref{sklar's Thm} and Lipschitz condition \eqref{Lipschitz}, we have
\begin{small}
\begin{align*}
I_{1}&=
\big|C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m}))-C_{R^{\ast}}(H^{\ast}_{i}(t_{1}),\ldots, H^{\ast}_{i}(t_{m}))\big|\\&\leq\sum_{i=1}^{m}\big| (F^{\ast}_{i}(t_{i})-H^{\ast}_{i}(t_{i})\big|.
\end{align*}
\end{small}
Note that, from the property of the Dirichlet process, for any $t_{i}\in\mathbb{R}$ and $\epsilon>0$, Chebyshev's inequality implies
\begin{align*}
Pr\left\lbrace \big| F^{\ast}_{i}(t_{i})-H^{\ast}_{i}(t_{i})\big|\geq\epsilon\right\rbrace
\leq\dfrac{H^{\ast}_{i}(t_{i})\left(1-H^{\ast}_{i}(t_{i})\right)}{(a+n+1)\epsilon^{2}}.
\end{align*}
Substituting $n=k^{2}+b$, for $k\in\mathbb{N}$ and $b\in\lbrace 0,1,\ldots\rbrace$, gives
\begin{align*}
Pr\left\lbrace \big| F^{\ast}_{i}(t_{i})-H^{\ast}_{i}(t_{i})\big|\geq\epsilon\right\rbrace
\leq\dfrac{1}{4\,k^{2}\epsilon^{2}}.
\end{align*}
Since $\sum_{k=1}^{\infty}k^{-2}$ converges, then $\sum_{k=1}^{\infty}Pr\left\lbrace \big| F^{\ast}_{i}(t_{i})-H^{\ast}_{i}(t_{i})\big|\geq\epsilon\right\rbrace<\infty.$ Hence, by the first Borel Cantelli lemma, $\big| F^{\ast}_{i}(t_{i})-H^{\ast}_{i}(t_{i})\big|\xrightarrow{a.s.}0$, as $k\rightarrow\infty$ or $n\rightarrow\infty$. This completes the proof.
\end{proof}
\section{Selecting $a$, $H$ and the Method of Estimation of $R^{\ast}$ in the BSPGC Approach}
The proposed method for modeling multivariate distributions depends on $a$, $H$ and $R^{\ast}$. Hence, it is necessary to look carefully at the impact of these parameters on the approach. For instance, from \eqref{pos base measure}, a large value of $a$ can increase the effect of the $m$-variate distribution $H$ instead of $F_{n}$ in the posterior-based model \eqref{pos-BSPGC}. The following lemma shows the effect of the value of $a$ on the model \eqref{pos-BSPGC}.
\begin{lemma}\label{effect a-pos}
Let $\mathbf{x}_{m\times n}$ be a sample from $m$-variate distribution function $F_{true}$ with unknown marginal cdf's $F_{1},\ldots,F_{m}$. Also, let $H$ be a known $m$-variate cdf with marginal cdf's $H_{1},\ldots,H_{m}$ and $H^{\ast}$ be the $m$-variate cdf, defined in \eqref{pos base measure}, with marginal cdf's $H^{\ast}_{1},\ldots,H^{\ast}_{m}$. Assume that $F^{\ast}_{i}\sim DP(a+n, H^{\ast}_{i})$, for $1\leq i\leq m$. Then, for any $\mathbf{t}\in \mathbb{R}^{m}$, $C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m}))\xrightarrow{a.s.}H(\mathbf{t})$ as $a\rightarrow\infty$.
\end{lemma}
\begin{proof}
The proof is similar to the proof of Lemma \ref{post as n to inf}. For this, assume that $a=k^{2}c$ for $k\in\lbrace 0,1,\ldots\rbrace$ and a fixed positive number $c$. For any fixed $n$, replace $F_{true}(\mathbf{t})$ and $a$, respectively by, $H(\mathbf{t})$ and $k^{2}c$ in the proof of Lemma \ref{post as n to inf}. Then the result follows.
\end{proof}\\

It follows from Lemma \ref{effect a-pos} that increasing the value of $a$ can lead to some errors. To avoid this issue, we propose to choose $a$ to be at most $0.5\,n$ as recommended in Al-Labadi and Zarepur (2017).

The choice of $H$ is also very significant and there are two main issues to reflect. The first one is the independence of the approach to the choice of $H$ (invariance property).   As pointed in Table \ref{priori-invariant}, the approach is invariance to the choice of any continuous multivariate distribution. The second issue is the compatibility between $H$ and the data. This is typically called \emph{prior-data conflict} (Evans and Moshonov, 2006; Al-Labadi and Evans, 2017,  Al-Labadi and Evans, 2018, Al-Labadi and Wang, 2019). As illustrated in Section 7.2, the existence of prior-data conflict yields to a failure of the approach and thus should be avoided. Since $E(F_{i})=H_{i}$, where $F_{i}\sim DP(a,H_{i})$ for $i=1,\ldots,m$, a reasonable choice of $H$ that ensures the avoidance of prior-data conflict  is the $m$-variate normal distribution $N_{m}(\overline{\mathbf{x}},S_{\mathbf{x}})$, where $\overline{\mathbf{x}}=1/n\sum_{i=1}^{n}\mathbf{x}_{i}$ and $S_{\mathbf{x}}=1/(n-1)\sum_{i=1}^{n}(\mathbf{x}_{i}-\overline{\mathbf{x}})(\mathbf{x}_{i}-\overline{\mathbf{x}})^{T}$.

To carry on the approach, it is essential to estimate the correlation matrix $R^{\ast}$. For this, we first generate a sample from the mixture distribution in \eqref{pos base measure}. Then, based on the generated sample, $R^{\ast}$  is estimated by one of the following three common procedures: the Gaussian correlation rank, the Kendal's $\tau$ or the Spearman's $\rho$. In Section 7, we performed a simulation study to compare the effect of these three methods on the quality of the approach. As a result,  we recommend using Kendal's $\tau$ correlation coefficients with $a=1$ in the proposed model.

\section{A MVN Test Based on the BSPGC Approach}
Let $\mathbf{x}_{m\times n}$ be a sample of size $n$ from an unknown $m$-variate  distribution $F_{true}$. The  problem to be addressed in this section is to test the hypothesis
\begin{equation}\label{test1}
\mathcal{H}_{0}:F_{true}\in \mathcal{F},
\end{equation}
\noindent where $\mathcal{F}=\{N_{m}(\boldsymbol{\mu}_{m},\Sigma_{m}):\boldsymbol{\mu}_{m}\in\mathbb{R}^{m},\det (\Sigma_{m}) >0\}$. Note that, whenever $\boldsymbol{\mu}_{m}$ and $\Sigma_{m}$ are unknown, they are to be estimated by the sample mean vector $\overline{\mathbf{x}}$ and sample covariance matrix
$S_{\mathbf{x}}$, respectively. Thus, for $\boldsymbol{\theta}_{\mathbf{x}}=(\overline{\mathbf{x}},S_{\mathbf{x}})$, $F_{\boldsymbol{\theta}_{\mathbf{x}}}=N_{m}(\overline{\mathbf{x}},S_{\mathbf{x}})$
is the best representative of the family $\mathcal{F}$ to compare with distribution $F_{true}$. Hence, testing \eqref{test1} is equivalent to test
\begin{equation}\label{test2}
\mathcal{H}_{0}:F_{true}=F_{\boldsymbol{\theta}_{\mathbf{x}}}.
\end{equation}
Now,  we continue as follows.
Let $H=F_{\boldsymbol{\theta}_{\mathbf{x}}}$ with marginal cdf's $H_{1}=F_{\boldsymbol{\theta}_{\mathbf{x}_{1}}},\ldots,H_{m}=F_{\boldsymbol{\theta}_{\mathbf{x}_{m}}}$. Here, for $i=1,\ldots,m$, $F_{\boldsymbol{\theta}_{\mathbf{x}_{i}}}$ is the cdf of the univariate normal distribution with mean $\overline{x}_{i}$ and variance $s^{2}_{i}$, where $\overline{x}_{i}$ and $s^{2}_{i}$ are the $i$-th element of $\overline{\mathbf{x}}$ and $i$-th diagonal element of $S_{\mathbf{x}}$, respectively. Assume that $F_{i}\sim DP(a, H_{i})$. For any $\mathbf{t}\in\mathbb{R}^{m}$, let
\begin{equation}\label{pri-BSPGC}
F(t_{1},\ldots,t_{m})=C_{R_{\mathbf{x}}}(F_{1}(t_{1}),\ldots, F_{m}(t_{m})),
\end{equation}
be the \emph{prior-based model}, where $R_{\mathbf{x}}=\left(r_{\mathbf{x},ij}\right)$ is the correlation matrix of the $m$-variate distribution $F_{\boldsymbol{\theta}_{\mathbf{x}}}$ and to be estimated as discussed in Section 4. Note that, as pointed out earlier, setting  $H_{i}=F_{\boldsymbol{\theta}_{\mathbf{x}_{i}}}$ ensures compatibility between the data and the prior which will certainly avoid prior-data conflict. More details about the effect of the prior-data conflict on the approach is clarified in Section 7, where it is revealed that  the existence of prior data conflict  leads to erroneous result of the test.

Recalling  the posterior-based model as defined in Section 3, to proceed  with the approach, the energy distance is used to compute the distance between this model and $F_{\boldsymbol{\theta}_{\mathbf{x}}}$ (posterior distance) and the distance between the prior-based model and $F_{\boldsymbol{\theta}_{\mathbf{x}}}$ (prior distance). The next lemma proposes a Bayesian counterpart of the distance \eqref{E-distance-2} as an appropriate tool to measure dissimilarities between the proposed models and the null distribution. This is considered a very convenient tool for assessing MVN in high dimensional problems ($m>n$).
\begin{lemma}\label{Bayesian-lemma}
Let $F_{true}$ be an $m$-variate distribution function with unknown marginal cdf's $F_{1},\ldots,F_{m}$ and $H$ be a known $m$-variate distribution function with marginal cdf's $H_{1},\ldots,H_{m}$. Assume that $F_{i}\sim DP(a,H_{i})$, and $F_{N_i}$'s are the approximation of Dirichlet process $F_{i}$'s, given by Ishwaran and Zarepour (2002). Consider the Energy distance between $F_{N}=C_{R}(F_{N_1}(t_{1}),\ldots,F_{N_m}(t_{m}))$ and $H$ as
\begin{equation}\label{Beyesian-E-distance}
d_{\mathcal{E},N,a}(F_{N},H)=2\sum_{i=1}^{N}J_{i,N}E_{H}||\mathbf{x}_{i}-\mathbf{Y}||-\sum_{i,j=1}^{N}J_{i,N}J_{j,N}||\mathbf{x}_{i}-\mathbf{x}_{j}||-E_{H}||\mathbf{Y}-\mathbf{Y}^{\prime}||,
\end{equation}
where $(J_{i,N})_{1\leq i\leq N}\sim Dirichlet(a/N,\ldots,a/N)$, $\mathbf{Y},\mathbf{Y}^{\prime}\overset{i.i.d.}{\sim}H$, $(\mathbf{x}_{1},\ldots,\mathbf{x}_{N})$ is an observed sample from $F_{N}$ and $R$ is the correlation matrix of $H$. Then, as $a\rightarrow\infty$
\begin{equation*}
E_{F_{N}}(d_{\mathcal{E},N,a}(F_{N},H))\xrightarrow{a.s.}d_{\mathcal{E},N}(F_{N},H),
\end{equation*}
where $d_{\mathcal{E},N}(F_{N},H)$ is defined in \eqref{E-distance-2} with $F=F_{n}$, $G=H$ and $n=N.$
\end{lemma}
\begin{proof}
From properties of Dirichlet distribution, $E_{F_{N}}(J_{i,N})=1/N$ and $E_{F_{N}}(J_{i,N}J_{j,N})$ $=a/((a+1)N^{2})$. Then
\begin{equation}\label{mean-E-dist}
E_{F_{N}}(d_{\mathcal{E},N,a}(F_{N},H))=\frac{2}{N}\sum_{i=1}^{N}E_{H}||\mathbf{x}_{i}-\mathbf{Y}||-\frac{a}{(a+1)N^{2}}\sum_{i,j=1}^{N}||\mathbf{x}_{i}-\mathbf{x}_{j}||-E_{H}||\mathbf{Y}-\mathbf{Y}^{\prime}||.
\end{equation}
The proof is immediately followed by letting $a\rightarrow\infty$ in \eqref{mean-E-dist}.
\end{proof}

The next lemma allows us to use the approximation of the Dirichlet process in the prior-based and posterior-based models for approximating the distribution of the posterior and the prior distances computed by \eqref{Beyesian-E-distance}.
\begin{lemma}
Let $F_{true}$ be an $m$-variate distribution function with unknown marginal cdf's $F_{1},\ldots,F_{m}$ and $H$ be a known $m$-variate distribution function with marginal cdf's $H_{1},\ldots,H_{m}$. Assume that $F_{i}\sim DP(a,H_{i})$ and $F_{Ni}$'s are the approximation of the Dirichlet process $F_{i}$'s, given in \eqref{approx of DP}. Then, for any $(t_{1},\ldots,t_{m})\in\mathbb{R}^{m}$, $C_{R}(F_{N1}(t_{1}),\ldots,$ $F_{Nm}(t_{m}))\xrightarrow{a.s.}C_{R}(F_{1}(t_{1}),\ldots, F_{m}(t_{m}))$, as $N\rightarrow\infty$, where $R$ is the correlation matrix of $H$.
\end{lemma}

\begin{proof}
From Lipschitz condition \eqref{Lipschitz}, we have
\begin{small}
\begin{align*}
\big|C_{R}(F_{N1}(x_{1}),\ldots, F_{Nm}(x_{m}))-C_{R}(F_{1}(x_{1}),\ldots, F_{m}(x_{m}))\big|\leq\sum_{i=1}^{m}\big| F_{N1}(x_{i})-F_{1}(x_{i})\big|.
\end{align*}
\end{small}
Since $F_{Ni}(x_{i})\xrightarrow{a.s.}F_{i}(x_{i})$, for $1\leq i\leq m$ (Zarepur and Al-Labadi, 2012), the result follows.
\end{proof}

The procedure is continued by considering $d_{\mathcal{E},N,a}(F_{N},F_{\boldsymbol{\theta}_{\mathbf{x}}})$ as the Energy distance between the prior-based model \eqref{pri-BSPGC} and the null distribution $F_{\boldsymbol{\theta}_{\mathbf{x}}}$ using formula \eqref{Beyesian-E-distance}. Similarly, consider $d_{\mathcal{E},N,a}(F^{\ast}_{N},F_{\boldsymbol{\theta}_{\mathbf{x}}})$ for the posterior-based model \eqref{pos-BSPGC} and the null distribution. Then, the relative belief ratio is used to compare the concentration of the posterior distribution  $d_{\mathcal{E},N,a}(F^{\ast}_{N},F_{\boldsymbol{\theta}_{\mathbf{x}}})$ to the prior distribution $d_{\mathcal{E},N,a}(F_{N},F_{\boldsymbol{\theta}_{\mathbf{x}}})$ about zero. As shown in the next lemma, if $\mathcal{H}_{0}$ is true, the distribution of the posterior distance should be more concentrated about 0 than the distribution of the prior distance; otherwise, the distribution of the prior distance should be more concentrated at 0 than the distribution of posterior distance. The comparison is made by computing the relative belief ratio with the interpretation as discussed
in$\ $Section 2.
\begin{lemma}\label{lem-2}
Let $\mathbf{x}_{m\times n}$ be a sample from $m$-variate distribution function $F_{true}$ with unknown marginal cdf's $F_{1},\ldots,F_{m}$. Assume that $\boldsymbol{\theta}_{\mathbf{x}}\xrightarrow{a.s.}\boldsymbol{\theta}_{0}$ and $F_{\boldsymbol{\theta}_{\mathbf{x}}}\xrightarrow{a.s.} F_{\boldsymbol{\theta}_{0}}$ as $n\rightarrow\infty$. Let $F^{\ast}_{i}\sim DP(a+n, H^{\ast}_{i})$, for $i=1,\ldots,m$. For any $\mathbf{t}=(t_{1},\ldots,t_{m})\in \mathbb{R}^{m}$ as $n\rightarrow\infty$

\begin{enumerate}[nolistsep,label=(\roman*),leftmargin=\parindent,align=left,labelwidth=\parindent,labelsep=0pt]
\item
$
\big|C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m}))-F_{\boldsymbol{\theta}_{\mathbf{x}}}(t_{1},\ldots,t_{m})\big|\xrightarrow{a.s.}0,
$ when $\mathcal{H}_{0}$ is true.
\item
$
\liminf\big|C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m}))-F_{\boldsymbol{\theta}_{\mathbf{x}}}(t_{1},\ldots,t_{m})\big|\displaystyle{\overset{a.s.}{>}} 0,
$ when $\mathcal{H}_{0}$ is not true,
\end{enumerate}
where $R^{\ast}$ is the correlation matrix of $H^{\ast}$, defined in \eqref{pos base measure}.
\end{lemma}

\begin{proof}
To prove (i), substitute $F_{true}$ in \eqref{triangular} by $F_{\boldsymbol{\theta}_{\mathbf{x}}}$. From \eqref{pos base measure}, for any $\mathbf{t}\in\mathbb{R}^{m}$, $H^{\ast}(\mathbf{t})\xrightarrow{a.s.} F_{true}(\mathbf{t})$ as $n\rightarrow\infty$. If $\mathcal{H}_{0}$ is true, then $F_{true}(\mathbf{t})=F_{\boldsymbol{\theta}_{0}}(\mathbf{t})$. Hence, the proof of (i) immediately follows from the proof of Lemma \ref{post as n to inf}. To prove (ii), Consider $I_{1}$ as in \eqref{triangular}. Applying the triangle inequality gives
\begin{align*}
\big|C_{R^{\ast}}(F^{\ast}_{1}(t_{1}),\ldots, F^{\ast}_{m}(t_{m}))-F_{\boldsymbol{\theta}_{\mathbf{x}}}(t_{1},\ldots,t_{m})\big|\geq \big|H^{\ast}(\mathbf{t})- F_{\boldsymbol{\theta}_{\mathbf{x}}}(\mathbf{t})\big|-I_{1}.
\end{align*}
Similar to the proof of Lemma \ref{post as n to inf}, $I_{1}\xrightarrow{a.s.}0$ and $\big|H^{\ast}(\mathbf{t})- F_{\boldsymbol{\theta}_{\mathbf{x}}}(\mathbf{t})\big|\xrightarrow{a.s.}\big|F_{true}(\mathbf{t})- F_{\boldsymbol{\theta}_{0}}(\mathbf{t})\big|$ as $n\rightarrow\infty$. Since $\mathcal{H}_{0}$ is not true, $\big|F_{true}(\mathbf{t})- F_{\boldsymbol{\theta}_{0}}(\mathbf{t})\big|\displaystyle{\overset{a.s.}{>}} 0$ which completes the proof of (ii).
\end{proof}

The effect of the value of $a$ on the posterior-based model was considered in Lemma \ref{effect a-pos}. It is also interesting to consider the effect of the value of $a$ on the proposed MVN test.
\begin{lemma}\label{effect a-pri}
Let $\mathbf{x}_{m\times n}$ be a sample from $m$-variate distribution function $F_{true}$ with unknown marginal cdf's $F_{1},\ldots,F_{m}$. Let $F_{\boldsymbol{\theta}_{\mathbf{x}}}$ be the cdf of $N_{m}(\overline{\mathbf{x}},S_{\mathbf{x}})$ with marginal cdf's $F_{\boldsymbol{\theta}_{\mathbf{x}_{1}}},\ldots,F_{\boldsymbol{\theta}_{\mathbf{x}_{m}}}$. If $F_{i}\sim DP(a,F_{\boldsymbol{\theta}_{\mathbf{x}_{i}}})$, for $i=1,\ldots,m$, then $C_{R_{\mathbf{x}}}(F_{1}(t_{1}),\ldots, F_{m}(t_{m}))$ $\xrightarrow{a.s.}F_{\boldsymbol{\theta}_{\mathbf{x}}}(\mathbf{t})$ as $a\rightarrow\infty$, where $R_{\mathbf{x}}$ is the correlation matrix of $F_{\boldsymbol{\theta}_{\mathbf{x}}}$.
\end{lemma}

\begin{proof}
The proof is similar to the proof of Lemma \ref{effect a-pos} and is omitted.
\end{proof}

Lemmas \ref{effect a-pos} and \ref{effect a-pri} show that for a too large value of $a$  (relative to the sample size) both the posterior-based and prior-based models are approaching to the null model $F_{\boldsymbol{\theta}_{\mathbf{x}}}$. Hence, the comparison between the posterior and prior distance to detect the normality can lead to an error in which we may accept $\mathcal{H}_{0}$ when it is not true and reject $\mathcal{H}_{0}$ when it is true. As recommended in Section 7, we should consider $a$ at most $0.5\,n$.

At the end of this section, it is worth pointing out that the proposed test can be extended to assess any family of multivariate distributions. For this, it is enough to consider a different family of multivariate distributions in the hypothesis \eqref{test1} and use its best representative distribution as $H$ in the methodology, which may be more challenging for some multivariate models.

\section{Main Steps for Testing the MVN}
The following computational algorithm summaries the main steps to test $\mathcal{H}_{0}$. This algorithm is viewed as a generalized version of Algorithm B of Al-Labadi and Evans (2018). Observe that, since closed forms of the densities of $D_{\mathcal{E}}=d_{\mathcal{E},N,a}(F_{N},F_{\boldsymbol{\theta}_{\mathbf{x}}})$ and  $D_{\mathcal{E}}|\mathbf{x}=d_{\mathcal{E},N,a}(F^{\ast}_{N},F_{\boldsymbol{\theta}_{\mathbf{x}}})$ are typically not available,   relative belief ratios need to be approximated via simulation.

\noindent\rule[0.1ex]{\linewidth}{1pt}\\[-.1cm]
\textbf{Algorithm 3} \vspace{-.05cm}{Relative belief algorithm based on the BSPGC approach for testing MVN}
\\[-.08cm]\noindent\rule[0.5ex]{\linewidth}{.4pt}
\begin{small}
\begin{enumerate}
\smallskip
\item Use Algorithm 2 to generate (approximately) marginal cdf $F_{i}$'s from $ DP(a, F_{\boldsymbol{\theta}_{\mathbf{x}_{i}}})$, for $1\leq i\leq m$.

\item Generate a sample of $N$ values from the $m$-variate distribution $F_{\boldsymbol{\theta}_{\mathbf{x}}}$ and estimate the correlation matrix $R_{\mathbf{x}}$, denoted by $\widehat{R}_{\mathbf{x}}$, as discussed in  Section 4.

\item Use the generated marginal cdf's $F_{i}$ and set $R=\widehat{R}_{\mathbf{x}}$ in Algorithm 1 to get a sample of $N$ values from prior-based model \eqref{pri-BSPGC}.

\item Use \eqref{Beyesian-E-distance} for the sample generated  in steps 3 to compute the prior distance $d_{\mathcal{E},N,a}(F_{N},F_{\boldsymbol{\theta}_{\mathbf{x}}})$.

\item Repeat steps (1)-(4) to obtain a sample of $r$ values from the
prior of $D_{\mathcal{E}}$.

\item Repeat steps (1)-(5) by replacing $a$ by $a+n$, $F_{i}$ by $F^{\ast}_{i}$,  $F_{\boldsymbol{\theta}_{\mathbf{x}_{i}}}$ by $H^{\ast}_{i}$, $F_{\boldsymbol{\theta}_{\mathbf{x}}}$  by $H^{\ast}$, $R_{\mathbf{x}}$ by $R^{\ast}$, $F_{N}$ by $F^{\ast}_{N}$, $D_{\mathcal{E}}$ by $D_{\mathcal{E}}|\mathbf{x}$ and prior by posterior. This yields to a sample of $r$ values from the
posterior of $D_{\mathcal{E}}|\mathbf{x}$.

\item Let $M$ be a positive number. Let $\hat{F}_{D_{\mathcal{E}}}$ denote the
empirical cdf of $D_{\mathcal{E}}$ based on the prior sample in step (5) and for $i=0,\ldots,M,$
let $\hat{d}_{i/M}$ be the estimate of $d_{i/M},$ the $(i/M)$-th prior
quantile of $D_{\mathcal{E}}.$ Here $\hat{d}_{0}=0$, and $\hat{d}_{1}$ is the largest value
of $d_{\mathcal{E}}$. Let $\hat{F}_{D_{\mathcal{E}}}(\cdot\,|\mathbf{x})$ denote the empirical cdf of $D_{\mathcal{E}}|\mathbf{x}$ based
on the posterior sample in step (10). For $d\in\lbrack\hat{d}_{i/M},\hat
{d}_{(i+1)/M})$, estimate $RB_{D_{\mathcal{E}}}(d\,|\,\mathbf{x})={\pi_{D_{\mathcal{E}}}(d|\mathbf{x})}/{\pi_{D_{\mathcal{E}}}(d)}$ by
\begin{equation}
\widehat{RB}_{D_{\mathcal{E}}}(d\,|\,\mathbf{x})=M\{\hat{F}_{D_{\mathcal{E}}}(\hat{d}_{(i+1)/M}\,|\,\mathbf{x})-\hat{F}%
_{D_{\mathcal{E}}}(\hat{d}_{i/M}\,|\,\mathbf{x})\}, \label{rbest}%
\end{equation}
the ratio of the estimates of the posterior and prior contents of $[\hat
{d}_{i/M},\hat{d}_{(i+1)/M}).$ Thus, we estimate $RB_{D_{\mathcal{E}}}(0\,|\,\mathbf{x})=\frac{\pi_{D_{\mathcal{E}}}(0|\mathbf{x})}{\pi_{D_{\mathcal{E}}}(0)}$
 by $\widehat{RB}_{D_{\mathcal{E}}}(0\,|\,\mathbf{x})=M\widehat{F}_{D_{\mathcal{E}}}(\hat{d}_{p_{0}}\,|\,\mathbf{x})$ where
$p_{0}=i_{0}/M$ and $i_{0}$ is chosen so that $i_{0}/M$ is not too small
(typically $i_{0}/M\approx0.05)$.\textbf{\smallskip}
\item Estimate the strength $DP_{D_{\mathcal{E}}}\big(RB_{D_{\mathcal{E}}}(d\,|\,\mathbf{x})\leq RB_{D_{\mathcal{E}}}%
(0\,|\,\mathbf{x})\,|\,\mathbf{x}\big)$ by the finite sum
\begin{equation}
\sum_{\{i\geq i_{0}:\widehat{RB}_{D}(\hat{d}_{i/M}\,|\,\mathbf{x})\leq\widehat{RB}%
_{D}(0\,|\,\mathbf{x})\}}\big(\hat{F}_{D_{\mathcal{E}}}(\hat{d}_{(i+1)/M}\,|\,\mathbf{x})-\hat{F}_{D_{\mathcal{E}}}(\hat
{d}_{i/M}\,|\,\mathbf{x})\big). \label{strest}%
\end{equation}
\end{enumerate}
\end{small}
\noindent\rule[0.5ex]{\linewidth}{.4pt}

For fixed $M,$ as $N\rightarrow\infty,r\rightarrow\infty,$
then $\hat{d}_{i/M}$ converges almost surely to $d_{i/M}$ and (\ref{rbest})
and (\ref{strest}) converge almost surely to $RB_{D_{\mathcal{E}}}(d\,|\,\mathbf{x})$ and
$DP_{D_{\mathcal{E}}}\big(RB_{D_{\mathcal{E}}}(d\,|\,\mathbf{x})\leq RB_{D_{\mathcal{E}}}(0\,|\,\mathbf{x})\,|\,\mathbf{x}\big)$, respectively.  The consistency of the proposed test is achieved by Proposition 6 of Al-Labadi and Evans (2018).  As a recommendation, one should try different values of $a$ to make sure the right conclusion has been obtained. However, we found out that setting $a=1$ gives adequate results. More details about implementing the approach is discussed in the following section.

\section{Simulation Studies}
This section is divided into two subsections. In the first subsection, the quality of the approach to model multivariate distributions is investigated, where different choices of $a$, $H$ and $R^{\ast}$ are considered. The evaluation technique relies on using the mean of the Energy distance $\overline{d}_{\mathcal{E},N}(F^{\ast},F_{true})$ based on $r$ replications. Note that, from Lemma \ref{Bayesian-lemma}, one may consider using  the package \textbf{energy} available in $\mathsf{R}$ to compute the distance.  We generated samples each of size $n=1000$ from a variety of bivariate distributions. The notations of the used distributions are listed (Table \ref{notation}) in Appendix A. In this study, we set $N=r=1000$ in Algorithm 3 with steps (1)-(6). Note that, for the methodology to work well, we expect $\overline{d}_{\mathcal{E},N}(F^{\ast},F_{true})$ to be close to zero. In the  second subsection, the proposed test is illustrated through several examples.

\subsection{Checking the Quality of the Posterior-based Model}
The performance of the posterior-based model (i.e. the quality of estimating the model) is illustrated by  considering the bivariate distributions given in Table \ref{performance} with some choices of $a$. The results are reported based on the Kendall's correlation coefficients. From Table \ref{performance}, the close values of $\overline{d}_{\mathcal{E},N}(F^{\ast},F_{true})$ to zero indicates to the good performance of the methodology to model bivariate distributions, particularly when  $a=1$. Note that, as mentioned in Lemma \ref{effect a-pos}, with increasing the value of $a$, the accuracy of the methodology will be decreased. For more illustration, part (a) of Figure \ref{box.a.norm} gives the boxplots of the energy distance between $F_{true}=N_{2}(\mathbf{0}_{2},A_{2})$ and its corresponding $F^{\ast}$ for $a=1,\,5,\,10.$ Boxplots of marginal distributions are also given in part (b) of Figure \ref{box.a.norm}. Also, the marginal densities of $N_{2}(\mathbf{0}_{2},A_{2})$ and its $F^{\ast}$ are given in Figure \ref{fig 1}. The bivariate scatter plots are shown below the diagonal, histograms on the diagonal and the Kendall correlation above the diagonal. Correlation ellipses and loess smooths (red lines) are also shown.
\setlength{\extrarowheight}{.1mm}
\begin{table}[ht]
\centering
\setlength{\tabcolsep}{2.1 mm}
\caption{The mean of the Energy distance between the true and the posterior-based model based on the Kendall's $\tau$ for $a=1,\,5,\,10$ and $n=1000$.}\label{performance}
\scalebox{0.86}{
\begin{tabular}{ccc||ccc}
\toprule
\bfseries True distribution & $a$ &\scalebox{1.2}{$\overline{d}_{\mathcal{E},N}(F^{\ast},F_{true})$}&\bfseries True distribution & $a$ &\scalebox{1.2}{$\overline{d}_{\mathcal{E},N}(F^{\ast},F_{true})$}\\\hline

&\scalebox{1.1}{1}& \scalebox{1.1}{0.00524} & & \scalebox{1.1}{1}&\scalebox{1.1}{0.00528} \\
\scalebox{1.1}{$N_{2}(\mathbf{0}_{2},A_{2})$} &\scalebox{1.1}{5}& \scalebox{1.1}{0.00559}  &\scalebox{1.1}{$(\mathcal{P}_{VII}(1,1,1))^{2}$} & \scalebox{1.1}{5} & \scalebox{1.1}{0.00551}\\
&\scalebox{1.1}{10}& \scalebox{1.1}{0.00581} & & \scalebox{1.1}{10} & \scalebox{1.1}{0.00573}\\
\hline
&\scalebox{1.1}{1}& \scalebox{1.1}{0.00517}  & & \scalebox{1.1}{1}&\scalebox{1.1}{0.00533}\\
\scalebox{1.1}{$t_{5}(\mathbf{0}_{2},I_{2})$} &\scalebox{1.1}{5}& \scalebox{1.1}{0.00538} &\scalebox{1.1}{$E(0.5)\otimes E(0.25)$} & \scalebox{1.1}{5} & \scalebox{1.1}{0.00575}  \\
&\scalebox{1.1}{10}& \scalebox{1.1}{0.00563} & & \scalebox{1.1}{10}&\scalebox{1.1}{0.00581} \\
\hline
&\scalebox{1.1}{1}& \scalebox{1.1}{0.00564} & & \scalebox{1.1}{1}& \scalebox{1.1}{0.00520} \\
\scalebox{1.1}{$LN_{2}(\mathbf{0}_{2},B_{2})$} &\scalebox{1.1}{5}& \scalebox{1.1}{0.00568} &\scalebox{1.1}{$B(1,2)\otimes B(2,1)$} & \scalebox{1.1}{5} & \scalebox{1.1}{0.00557} \\
&\scalebox{1.1}{10}& \scalebox{1.1}{0.00597} & & \scalebox{1.1}{10} & \scalebox{1.1}{0.00563}\\
\bottomrule
\end{tabular}
}
\end{table}

\begin{figure}[ht]
 \centering
    \subfloat[]{{\includegraphics[width=5.2cm]{box-a.pdf} }}%
    \qquad
    \subfloat[]{{\includegraphics[width=5.2cm]{box-group-a-normal.pdf} }}
    \caption{a: Boxplots of the Energy distance between $F_{true}=N_{2}(\mathbf{0}_{2},A_{2})$ and $F^{\ast}$ for $a=1,\,5,\,10$ and $n=1000$. b: Boxplots of the marginal distributions.}%
    \label{box.a.norm}%
\end{figure}
\begin{figure}[ht]
    \centering
    \subfloat[$F_{true}$: True Distribution.]{{\includegraphics[width=5.2cm]{true-norm-dot.pdf} }}%
    \qquad
    \subfloat[$F^{\ast}$: Posterior-based model.]{{\includegraphics[width=5.2cm]{post-norm-dot1.pdf} }}
    \caption{Marginal densities of $F_{true}=N_{2}(\mathbf{0}_{2},A_{2})$ and its posterior-based model with $n=1000$ based on the Kendall's $\tau$ for $a=1$.}%
    \label{fig 1}%
\end{figure}

Next, we inspect the effect of choosing different correlation coefficients such as the Gaussian rank, the Kendall's $\tau$ and the Spearman's $\rho$ on the posterior-based model. Consider $(\mathcal{P}_{VII}(1,1,1))^{2}$ and $N_{2}(\mathbf{0}_{2},A_{2})$ as two true distributions (consult  Table \ref{notation} in Appendix A for the notations). Table \ref{correlation} reports the results for $a=1$. Note that, the $\mathsf{R}$ package \textbf{rococo} is used to estimate $R^{\ast}$ based on the Gaussian rank correlation coefficients. It follows from Table \ref{correlation} that the performance of the methodology is approximately the same for different correlation coefficients.

\setlength{\extrarowheight}{.1mm}
\begin{table}[ht]
\centering
\setlength{\tabcolsep}{3 mm}
\caption{$\overline{d}_{\mathcal{E},N}(F^{\ast},F_{true})$ based on the Gaussian rank, Kendall, and Spearman correlation coefficients for $n=1000$ and $a=1$.}\label{correlation}
\scalebox{0.86}{
\begin{tabular}{c|ccc}
\toprule
\bfseries True distribution &\bfseries Gaussian rank &\bfseries Kendall's $\tau$ &\bfseries Spearman's $\rho$  \\
\hline
\scalebox{1.1}{$(\mathcal{P}_{VII}(1,1,1))^{2}$} &\scalebox{1.1}{0.00557} &\scalebox{1.1}{0.00528} &\scalebox{1.1}{0.00567}  \\
\scalebox{1.1}{$N_{2}(\mathbf{0}_{2},A_{2})$} &\scalebox{1.1}{0.00542} &\scalebox{1.1}{0.00524} &\scalebox{1.1}{0.00556} \\
\bottomrule
\end{tabular}
}
\end{table}
\subsection{Checking MVN Based on the BSPGC Approach}
The proposed normality test is illustrated through some interesting examples discussed in Henze and Visagie (2019). Note that, $NMIX1$ is a skewed heavy-tailed and $NMIX2$ is a symetric heavy-tailed distribution. Also, $(\mathcal{P}_{VII}(1,1,r))^{2}$, for $r\geq 10$ is a symetric distribution and has very similar behavior with a bivariate normal distribution. For a given sample of size  $n=50$, generated from distributions in Table \ref{T1}, the bivariate normality assumption is checked. For all cases, we set $N=r=1000$ and $M=20$ in Algorithm 3. To study the sensitivity of the approach, various values of $a$ are considered. The results of the proposed test are reported in Table \ref{T1}. The results are also compared to the Energy (E)-test (Sz\'{e}kely and Rizzo, 2013). Reminding that we want $RB>1$ and the strength close to 1 when $\mathcal{H}%
_{0}$ is true and $RB<1$ and the strength close to 0 when $\mathcal{H}_{0}$ is
false,  it is seen from  Table \ref{T2} that the proposed test has an excellent performance to accept or reject the bivariate normality assumption. The type I error and the power of the test are also reported in Table \ref{T2}. They show that the proposed test is powerful in both accepting and rejecting $\mathcal{H}_{0}$.

The next example uses a real data set.

\noindent\textbf{Real data example} (\emph{Swiss Heads}): In this example, we consider the data of six readings on the dimensions of the heads of 200 twenty year old Swiss soldiers given by Flury and Riedwyl (1988). The variables are minimal frontal breadth, breadth of angulus mandibulae, true facial height, length from glabella to apex nasi, length from tragion to nasion, and length from tragion to gnathion. The problem is to assess the six-variate normality assumption for this data set. The E-test's p-value is $2.2\times10^{-16}$, which shows strong evidence to reject the six-variate normality assumption. The proposed test presents $RB=0$ and strength$=0$ based on the Kendall's $\tau$ and $a=1$ which follows the methodology also presents strong evidence to reject the six-variate normality assumption.

We end this subsection by investigating   the effect of the prior-data conflict on the approach. This is in fact highlights the effect of the choice of $H$ in the prior-based model. For this, consider the results of MVN test when $F_{true}=(Exp(0.5))^{2}$ for different choices of $H$ in Table \ref{prior-data}. Clearly, when $H=F_{\boldsymbol{\theta}_{\mathbf{x}}}$ the results are correct; otherwise,  they  are incorrect.  Another concern is to check the effect of the double use of the data by considering $H$ as $F_{\boldsymbol{\theta}_{\mathbf{x}}}$ in the prior distance. Particularly, Table \ref{priori-invariant} gives the mean of the prior distance $\overline{d}_{\mathcal{E},N,a}(F_{N},H)$ for various choices of $H$. It is obvious from this table that the prior distance is invariant with respect to the choice of $H$
\setlength{\extrarowheight}{0.1 mm}
\begin{table}[h!]
\centering
\setlength{\tabcolsep}{.9 mm}
\caption{Relative belief ratios and strength (str)-s for testing the bivariate normality assumption
with various alternatives and choices of $a$ based on the Kendall's $\tau$ with $n=50$.}\label{T1}
\scalebox{0.86}{
\begin{tabular}{ccccc||ccccc}
\toprule
{\parbox{22mm}{\center{\bfseries{Alternative\hspace{4mm} distribution}}}} &\bfseries $a$& $RB$&\bfseries  $ str$&{\parbox{22mm}{\center{\bfseries{E-test's\hspace{4mm} p-value}}}}&{\parbox{22mm}{\center{\bfseries{Alternative\hspace{4mm} distribution}}}} &\bfseries $a$& $RB$&\bfseries  $ str$&{\parbox{22mm}{\center{\bfseries{E-test's\hspace{4mm} p-value}}}}\\[.43 cm] \hline
 &\scalebox{1.1}{1} & \scalebox{1.1}{3.54} & \scalebox{1.1}{0.823} & &&\scalebox{1.1}{1} & \scalebox{1.1}{0.04} & \scalebox{1.1}{0.005} & \\[-.07cm]
\scalebox{1.1}{$N_{2}(\mathbf{0}_{2},I_{2})$}&\scalebox{1.1}{5}& \scalebox{1.1}{3.26} & \scalebox{1.1}{0.832} &\scalebox{1.1}{0.8794}&\scalebox{1.1}{$E(0.5)\otimes E(0.25)$}&\scalebox{1.1}{5} & \scalebox{1.1}{0.00} & \scalebox{1.1}{0.000} &\scalebox{1.1}{$2.2\times10^{-16}$}\\[-.07cm]
&\scalebox{1.1}{10}  & \scalebox{1.1}{2.48} & \scalebox{1.1}{0.997}  & &&\scalebox{1.1}{10}  & \scalebox{1.1}{0.00} & \scalebox{1.1}{0.000} & \\[-.01cm]\hline

 & \scalebox{1.1}{1}& \scalebox{1.1}{3.76} & \scalebox{1.1}{0.999}& &  &\scalebox{1.1}{1} & \scalebox{1.1}{0.80} & \scalebox{1.1}{0.110} & \\[-.07cm]
\scalebox{1.1}{$N_{2}(\mathbf{0}_{2},A_{2})$}&\scalebox{1.1}{5} &  \scalebox{1.1}{2.82}  & \scalebox{1.1}{0.859} & \scalebox{1.1}{0.8442}&\scalebox{1.1}{$\mathcal{S}^{2}(LN(0,0.25))$}&\scalebox{1.1}{5} & \scalebox{1.1}{0.62} & \scalebox{1.1}{0.151} &\scalebox{1.1}{$2.2\times10^{-16}$}\\[-.07cm]
&\scalebox{1.1}{10} & \scalebox{1.1}{2.22}      & \scalebox{1.1}{0.884} & &&\scalebox{1.1}{10}  & \scalebox{1.1}{0.60} & \scalebox{1.1}{0.240}  & \\[-.01cm]\hline

&\scalebox{1.1}{1} & \scalebox{1.1}{0.18} & \scalebox{1.1}{0.017} && &\scalebox{1.1}{1}  & \scalebox{1.1}{2.92} & \scalebox{1.1}{0.854} & \\[-.07cm]
\scalebox{1.1}{$LN_{2}(\mathbf{0}_{2},B_{2})$}&\scalebox{1.1}{5} & \scalebox{1.1}{0.18} & \scalebox{1.1}{0.033} &\scalebox{1.1}{$2.2\times10^{-16}$}&\scalebox{1.1}{$(\mathcal{P}_{VII}(1,1,10))^{2}$}&\scalebox{1.1}{5} & \scalebox{1.1}{2.74} & \scalebox{1.1}{0.999} &\scalebox{1.1}{0.7035}\\[-.07cm]
&\scalebox{1.1}{10}  & \scalebox{1.1}{0.06} & \scalebox{1.1}{0.007}  && &\scalebox{1.1}{10} & \scalebox{1.1}{2.26} & \scalebox{1.1}{0.878}  & \\[-.01cm]\hline

&\scalebox{1.1}{1} & \scalebox{1.1}{0.10} & \scalebox{1.1}{0.000} & &&\scalebox{1.1}{1} & \scalebox{1.1}{0.34} & \scalebox{1.1}{0.017} & \\[-.07cm]
\scalebox{1.1}{$NMIX1$}&\scalebox{1.1}{5} & \scalebox{1.1}{0.02} & \scalebox{1.1}{0.002} &\scalebox{1.1}{0.0050}&\scalebox{1.1}{$(\chi^{2}_{5})^{2}$}&\scalebox{1.1}{5} & \scalebox{1.1}{0.22} & \scalebox{1.1}{0.026} &\scalebox{1.1}{0.03518}\\[-.07cm]
&\scalebox{1.1}{10}  & \scalebox{1.1}{0.00} & \scalebox{1.1}{0.000}  & &&\scalebox{1.1}{10} & \scalebox{1.1}{0.04} & \scalebox{1.1}{0.002}  & \\[-.01cm]\hline

&\scalebox{1.1}{1} & \scalebox{1.1}{3.64} & \scalebox{1.1}{1.000} && &\scalebox{1.1}{1} & \scalebox{1.1}{0.58} & \scalebox{1.1}{0.099} & \\[-.07cm]
\scalebox{1.1}{$NMIX2$}&\scalebox{1.1}{5} & \scalebox{1.1}{3.32} & \scalebox{1.1}{0.834} &\scalebox{1.1}{0.8744}&\scalebox{1.1}{$N(0,1)\otimes \chi^{2}_{5}$}&\scalebox{1.1}{5} & \scalebox{1.1}{0.40} & \scalebox{1.1}{0.070} &\scalebox{1.1}{0.4020}\\[-.07cm]
&\scalebox{1.1}{10}  & \scalebox{1.1}{2.34} & \scalebox{1.1}{0.874} &&&\scalebox{1.1}{10} & \scalebox{1.1}{0.30} & \scalebox{1.1}{0.095}  & \\[-.01cm]\hline

&\scalebox{1.1}{1} & \scalebox{1.1}{0.80} & \scalebox{1.1}{0.210} && &\scalebox{1.1}{1} & \scalebox{1.1}{0.63} & \scalebox{1.1}{0.160} & \\[-.07cm]
\scalebox{1.1}{$t_{5}(\mathbf{0}_{2},I_{2})$}&\scalebox{1.1}{5} & \scalebox{1.1}{0.60} & \scalebox{1.1}{0.181} &\scalebox{1.1}{0.0452}&\scalebox{1.1}{$N(0,1)\otimes t_{3}$}&\scalebox{1.1}{5} & \scalebox{1.1}{0.44} & \scalebox{1.1}{0.092} &\scalebox{1.1}{0.0502}\\[-.07cm]
&\scalebox{1.1}{10} & \scalebox{1.1}{0.20} & \scalebox{1.1}{0.010}  & &&\scalebox{1.1}{10} & \scalebox{1.1}{0.40} & \scalebox{1.1}{0.020}  & \\[-.01cm]\bottomrule

\end{tabular}
}
\end{table}

\setlength{\extrarowheight}{.1mm}
\begin{table}[ht]
\centering
\setlength{\tabcolsep}{7.5 mm}
\caption{The proportion of rejecting (POR) $\mathcal{H}_{0}$ out of 1000 replications based on the Kendall's $\tau$ for $a=1$ and sample of size $n=50$.}\label{T2}
\scalebox{0.86}{
\begin{tabular}{cc||cc}
\toprule\\[-.4cm]
\bfseries Distribution & \bfseries POR $\mathcal{H}_{0}$&\bfseries Distribution&\bfseries POR $\mathcal{H}_{0}$\\\hline
\scalebox{1.2}{$N_{2}(\mathbf{0}_{2},I_{2})$}&\scalebox{1.2}{0.057}$^{\dagger}$& \scalebox{1.2}{$NMIX1$} &\scalebox{1.2}{0.794}$^{\ddagger}$ \\

\scalebox{1.2}{$N_{2}(\mathbf{0}_{2},A_{2})$}&\scalebox{1.2}{0.070$^{\dagger}$}&  \scalebox{1.2}{$NMIX2$}  & \scalebox{1.2}{0.077$^{\ddagger}$} \\

\scalebox{1.2}{$LN_{2}(\mathbf{0}_{2},B_{2})$}&\scalebox{1.2}{0.801$^{\ddagger}$}&  \scalebox{1.2}{$(\mathcal{P}_{VII}(1,1,10))^{2}$}  & \scalebox{1.2}{0.106$^{\ddagger}$} \\

\scalebox{1.2}{$\mathcal{S}^{2}(LN(0,0.25))$}&\scalebox{1.2}{0.798$^{\ddagger}$}& \scalebox{1.2}{$(\chi^{2}_{5})^{2}$} & \scalebox{1.2}{0.782$^{\ddagger}$} \\

\scalebox{1.2}{$t_{5}(\mathbf{0}_{2},I_{2})$}&\scalebox{1.2}{0.499$^{\ddagger}$}&  \scalebox{1.2}{$N(0,1)\otimes \chi^{2}_{5}$}  & \scalebox{1.2}{0.657$^{\ddagger}$} \\

\scalebox{1.2}{$E(0.5)\otimes E(0.25)$}&\scalebox{1.2}{0.999$^{\ddagger}$}&  \scalebox{1.2}{$N(0,1)\otimes t_{3}$}  & \scalebox{1.2}{0.551$^{\ddagger}$} \\\bottomrule

\end{tabular}
}
\centering
\begin{tablenotes}
      \small
      \item \fontsize{8.6}{8}\selectfont{$^{\dagger}$ Type I error of the proposed test.}
      \item $^{\ddagger}$ Power of the proposed test.
    \end{tablenotes}
\end{table}
\setlength{\extrarowheight}{0 mm}
\begin{table}[H]
\centering
\setlength{\tabcolsep}{5 mm}
\caption{RB(Strength) of a sample of size 50 generated from $(E(0.5))^{2}$ when there is prior-data conflict (a tiny overlap between the effective support regions).}\label{prior-data}
\scalebox{0.86}{
\begin{tabular}{c|ccc}
\toprule
$H$ &$a=1$ &$a=5$ &$a=10$  \\
\hline
\scalebox{1.1}{$F_{\boldsymbol{\theta}_{\mathbf{x}}}$} &\scalebox{1.1}{0.080(0.004)} &\scalebox{1.1}{0.040(0.006)} &\scalebox{1.1}{0.00(0.000)}  \\
\scalebox{1.1}{$N_{2}(\mathbf{x},I_{2})$} &\scalebox{1.1}{18.46(1.000)} &\scalebox{1.1}{9.180(0.989)} &\scalebox{1.1}{9.08(0.969)} \\
\scalebox{1.1}{$N_{2}(\mathbf{3}_{2},S_{\mathbf{x}})$} &\scalebox{1.1}{19.00(1.000)} &\scalebox{1.1}{10.71(1.000)} &\scalebox{1.1}{8.32(0.582)} \\
\bottomrule
\end{tabular}
}
\end{table}
\setlength{\extrarowheight}{.5mm}
\begin{table}[ht]
\centering
\begin{small}
\setlength{\tabcolsep}{5 mm}
\caption{The mean of the prior distance $\overline{d}_{\mathcal{E},N,a}(F_{N},H)$ for various choices of $H$ when $\mathbf{x}_{2\times 50}\sim (E(0.5))^{2}$ based on the Kendall's $\tau$ with $N=1000$ and $a=1$.}\label{priori-invariant}
\scalebox{.86}{
\begin{tabular}{cc||cc}
\toprule
\scalebox{1.2}{$H$} &\scalebox{1.2}{$\overline{d}_{\mathcal{E},N,a}(F,H)$} & \scalebox{1.2}{$H$} & \scalebox{1.2}{$\overline{d}_{\mathcal{E},N,a}(F,H)$} \\
\hline
\scalebox{1.2}{$F_{\boldsymbol{\theta}_{\mathbf{x}}}$} & \scalebox{1.2}{0.00672} & \scalebox{1.2}{$N_{2}(\mathbf{3}_{2},A_{2})$} & \scalebox{1.2}{0.00633}  \\
\scalebox{1.2}{$N_{2}(\mathbf{0}_{2},I_{2})$} & \scalebox{1.2}{0.00661} & \scalebox{1.2}{$NMIX2$} & \scalebox{1.2}{0.00658} \\
\bottomrule
\end{tabular}
}
\end{small}
\end{table}

\section{Concluding Remarks}
A BSPGC approach and its application to the MVN test have been suggested. In this procedure, a Gaussian copula model has been utilized to induce the dependence structure of the underlying multivariate distribution $F_{true}$. The Dirichlet process then has been constructed on the unknown margins of $F_{true}$ to define the prior-based and posterior-based models, respectively. The test has been developed by using the relative belief ratio for comparing the concentration of the distribution of the distance between the posterior-based model and the null distribution versus the concentration of the distribution of the distance between the prior-based model and the null distribution at zero. The Energy distance has been applied to compute distances as an appropriate tool especially in high dimensional problems. The methodology has been examined by a simulation study to clarify its excellent performance. Finally, application of the test including a real data example has been presented. A main advantage of the procedure is that it takes into account the dependence structure of the data in the MVN test. The extension of the procedure to different areas of the multivariate data analysis by considering various families of copula will be a part of a future research work.


\section{Introduction}
Medical image registration
consists of finding a spatial transformation that aligns two or more images.
The intrinsic ill-posedness of registration can lead to anatomically implausible transformations
associated with unrealistic volumetric distortion of the anatomy.
For certain anatomical regions, such as the myocardium, physiologically plausible image registration requires volumetric preservation which corresponds to so-called \emph{incompressible} registration~\cite{Bistoquet2008,Mansi2011,Tobon-Gomez2013}.

Incompressible medical image registration between a pair of images $I_1, I_2$
can be defined as a constrained optimisation problem:
\begin{equation}
    \label{eq:reg_incompressible_continuous}

        \argmin_{\Phi \in \mathcal{D}\left(\Omega\right)}\,\,
        \mathcal{L}(I_1, I_2, \Phi) + R(\Phi) \quad
        \text{s.t.} \quad
        \left[
        \forall x \in \mathcal{M},\quad \det\left(J_{\Phi}(x)\right) = 1
        \right]

\end{equation}
where $\Omega$ is the image domain, $\mathcal{M} \subset \Omega$ is the incompressible region, $\mathcal{L}$ an image dissimilarity measure, $R$ a regularisation term, $\mathcal{D}\left(\Omega\right)$ the group of diffeomorphic transformations from $\Omega$ onto itself, and $J_{\Phi}$ the Jacobian matrix of the transformation $\Phi$.

In practice, to solve \eqref{eq:reg_incompressible_continuous}, the transformation must be parameterised by a finite number of parameters, and the constraint $\forall x \in \mathcal{M},\, \det\left(J_{\Phi}(x)\right) = 1$ must be reduced to a finite number of equality constraints.
There are two approaches to discretise the constraint:
1) relaxing the constraint into an additive soft constraint~\cite{Aganj2015,DeCraene2012,Heyde2016,Rohlfing2003}.
2) using a specific parameterisation for the transformation~\cite{Bistoquet2008,Mang2018,Mansi2011}.
Soft constraints introduce hyperparameters that are difficult to tune reliably and cannot guarantee an incompressible transformation.
In \cite{Bistoquet2008}, the transformation is parameterised by a divergence-free displacement. Yet, this approach only provides a first order approximation of an incompressible deformation, still requiring a soft constraint.
Recently \cite{Mang2018,Mansi2011} proposed
to parameterise the transformation by a divergence-free stationary velocity field (SVF).
The transformation is obtained via the Lie exponential mapping $\exp: v \mapsto \Phi $ that maps any divergence-free SVF $v$ into an incompressible transformation $\Phi$~\cite{Mang2018,Mansi2011}.
This reduces the non-convex constraint in \eqref{eq:reg_incompressible_continuous} to a linear constraint in the continuous domain.
However, in~\cite{Mang2018,Mansi2011} the SVF is parameterised as linear B-splines, and the linear constraint is discretised by imposing the constraint only on the points of the deformation grid.
As a result, guarantees to obtain a continuous incompressible transformation do not hold anymore.
To mitigate this issue,~\cite{Mang2018} proposed to work with images of higher resolution with a finer grid for the linear B-spline leading to the need for distributed super-computing.
Moreover,~\cite{Mang2018,Mansi2011} methods are limited to using sum of squared differences (SSD) as image similarity metric, which limits their applicability to images with similar intensity distributions.

In this paper, we propose a constrained optimisation framework for incompressible diffeomorphic registration that allows to use any smooth image similarity and regularisation penalty.
As an efficient means of discrete SVF parameterisation for this problem, we introduce multivariate divergence-conforming B-splines that have recently raised interest in computational physics~\cite{Evans2013a}.
We demonstrate that their properties can be exploited to impose bounds on the divergence of the SVF over the entire continuous space using sparse linear constraints on its finite parameters.
Our general problem formulation allows us to solve incompressible registration using any state-of-the-art optimiser for constrained non-convex optimisation (e.g. \texttt{IPOPT}~\cite{Wachter2006}).
For evaluation, we initially apply our method for multi-modal incompressible registration of synthetically deformed brains.
We then compare our approach against the state-of-the-art results published for the STACOM'11 myocardial tracking challenge dataset~\cite{Tobon-Gomez2013} and achieve similar results while better retaining the incompressibility.

\newpage
\section{Method}
\subsubsection{Divergence-conforming B-splines}
General B-splines are very popular to parameterise deformations and velocity fields
over a continuous spatial domain $\Omega$
with a finite number of parameters
$\phi_i \in \mathds{R}^{3}$:
\begin{equation}
    \label{eq:bspline3d}
     \forall \left(x, y, z\right) \in \Omega,\quad
     v\left(x,y,z\right) = \sum_{i} B^k_{i,X}(x) B^k_{i,Y}(y) B^k_{i,Z}(z)\,\phi_{i},
\end{equation}
where
$(x_i, y_i, z_i)$ are the knots of a regular grid of spacing $\delta x \times \delta y \times \delta z$ and
$B^k_{i,U}$ are the
1D B-spline basis functions of order $k$ in the direction $U\in \{X, Y, Z\}$ (see supplemental material for more details).
Popular choices of the order are $k=1$ (linear B-splines)~\cite{Mang2018,Mansi2011} and $k=3$ (cubic B-splines)~\cite{Tobon-Gomez2013}.
A fundamental property of 1D B-spline basis functions of the variable $u$, on a regular grid
of spacing $\delta u$ with $n$ knots is that:
\begin{equation}
\label{eq:deriv_bsplines}

    \forall i,
    \quad \dv{B^k_{i,U}}{u} = \frac{B^{k-1}_{i,U} - B^{k-1}_{i+1,U}}{\delta u}
\end{equation}
This implies that the derivative of a 1D B-spline on a regular grid is also a 1D B-spline on the same grid, albeit of lower order.
However, this property does not extend to 3D B-splines using definition \eqref{eq:bspline3d}.
Especially, the divergence of those 3D B-splines are not B-splines,
because of the mixed order appearing with first partial derivatives.
This limitation makes it more difficult to relate properties of the velocity field $v$ to the values of the parameters
$\phi_i$.

To overcome this limitation, we propose to use a 3D divergence-conforming B-splines~\cite{Evans2013a} to parameterise $v$.
The orders of the 3D B-splines basis of each component are chosen so that the divergence of $v$ is, in a continuous sense, exactly a 1D B-spline of order $k$.
Using the same notations as in \eqref{eq:bspline3d}, and
using $\phi_i = (\phi^X_i,\phi^Y_i,\phi^Z_i)$,
we have:
\begin{equation}
    \forall (x,y,z) \in \Omega,\quad
    v(x,y,z) =
    \left(
    \begin{array}{c}
         \sum_{i} B^{k+1}_{i,X}(x) B^k_{i,Y}(y) B^k_{i,Z}(z)\,\phi^X_{i}  \\
         \sum_{i} B^k_{i,X}(x) B^{k+1}_{i,Y}(y) B^k_{i,Z}(z)\,\phi^Y_{i} \\
         \sum_{i} B^k_{i,X}(x) B^k_{i,Y}(y) B^{k+1}_{i,Z}(z)\,\phi^Z_{i}
    \end{array}
    \right)
\end{equation}
Using \eqref{eq:deriv_bsplines}, we obtain the continuous divergence $\div{v}$ for all $(x,y,z) \in \Omega$:
\begin{equation}
\label{eq:div_divergence-conforming}
    \begin{split}
        \div{v}(x,y,z) =& \sum_{i} B^{k}_{i,X}(x) B^k_{i,Y}(y) B^k_{i,Z}(z)\, \psi_i \\
        \text{s.t.}\,\, \forall i,\quad \psi_i =& \,
    \frac{\phi^X_{i} - \phi^X_{i-1}}{\delta x} +
    \frac{\phi^Y_{i} - \phi^Y_{i-1}}{\delta y} +
    \frac{\phi^Z_{i} - \phi^Z_{i-1}}{\delta z}
    \end{split}
\end{equation}
As a consequence, the following lemma states that $\div{v}$ is uniformly bounded at any point of a continuous subregion $\mathcal{M} \subset \Omega$ provided that a finite number of linear constraints are satisfied by the coefficients
$(\phi^X_i,\phi^Y_i,\phi^Z_i)$.

\begin{lemma}
\label{lemma:1}
Let $k\geq 2$, and $\mathcal{M}$ be a non-empty subset of $\Omega$. Let $\epsilon \geq 0$
and let $\mathbf{J}_{\mathcal{M}}=\{i
\,\,|\,\,
\left(
\supp B^{k}_{i,X} \times \supp B^{k}_{i,Y} \times \supp B^{k}_{i,Z}
\right)
\cap \mathcal{M} \neq \emptyset
\}$.
If
\[
    \forall i \in \mathbf{J}_{\mathcal{M}},\quad -\epsilon \leq
    \frac{\phi^X_{i} - \phi^X_{i-1}}{\delta x} +
    \frac{\phi^Y_{i} - \phi^Y_{i-1}}{\delta y} +
    \frac{\phi^Z_{i} - \phi^Z_{i-1}}{\delta z}
    \leq \epsilon,
\]
then at any point $m \in \mathcal{M}$, it holds that
$\,\abs{\div{v}(m)} \leq \epsilon$.
\end{lemma}

\begin{proof}
The proof follows from \eqref{eq:div_divergence-conforming} and the fact that the value of a B-spline at any point $m$ is bounded by the values of the B-spline coefficients associated to the knots that are close to $m$ (see supplementary material for more details).
\end{proof}

\subsubsection{Optimisation formulation and implementation}
In this section, we formulate the optimisation problem for diffeomorphic registration with the proposed incompressibility constraint on a subregion $\mathcal{M}$.
Using Lemma~\ref{lemma:1} and previous notations, our (symmetric) optimisation formulation of \eqref{eq:reg_incompressible_continuous} is:
\begin{equation}
    \label{eq:optim_formulation}
    \begin{split}
        \argmin_{\Theta=\{(\phi^X_i,\phi^Y_i,\phi^Z_i)\}_{i}} \,\, & \mathcal{L}(I_1 \circ \widetilde{\exp}(v(\Theta)), I_2) + \mathcal{L}(I_1, I_2 \circ \widetilde{\exp}(-v(\Theta))) + R(v(\Theta))\\
        \text{s.t. } \quad &
        \forall i \in \mathbf{J}_{\mathcal{M}},\,\,
            \frac{\phi^X_{i} - \phi^X_{i-1}}{\delta x} +
            \frac{\phi^Y_{i} - \phi^Y_{i-1}}{\delta y} +
            \frac{\phi^Z_{i} - \phi^Z_{i-1}}{\delta z}
            = 0\\
    \end{split}
\end{equation}
where
$\widetilde{\exp}$ is an approximation
of the Lie exponential.
Thus we obtain a constrained optimisation formulation that guarantees the SVF to be \emph{exactly} divergence-free over the entire continuous subregion $\mathcal{M}$ and that can be solved with efficient state-of-the-art optimisers.
Using state-of-the-art solvers like \texttt{IPOPT}, that uses a primal-dual interior-point filter line-search method, an approximated solution of \eqref{eq:optim_formulation} that satisfies the constraints up to machine precision can be obtained.

\subsubsection{Lie exponential approximation and incompressibility}
As the Lie exponential has to be approximated in practice, having a divergence-free velocity field does not strictly guarantee that the resulting deformation will be incompressible.
We now quantify this approximation in our framework with respect to the time step $\tau$ for the Euler method.
Let $v$ denote a (SVF) solution of \eqref{eq:optim_formulation} that fulfills the divergence-free constraint up to machine precision $\epsilon_{mach}$. For any point $m \in \mathcal{M}$, the Jacobian of the first step of the Euler method, $I + \tau v$, fulfills:
\begin{equation}
    \label{eq:jac_one_step_error}
    \det\left(J_{I + \tau v}(m)\right) = 1 + \tau \div{v}(m) + \mathcal{O}\left(\tau^2\right)
    \quad \text{where} \quad \abs{\div{v}(m)} \leq \epsilon_{mach}.
\end{equation}
Then, by composing $1/\tau$ times, if the point $m$ remains inside $\mathcal{M}$ during the Euler integration, the Jacobian of $\widetilde{\exp}(v)$ satisfies (see supplementary material):
\begin{equation}
    \label{eq:jac_error}
    \det\left(J_{\widetilde{exp}(v)}(m)\right) = 1 + \mathcal{O}\left(\tau + \epsilon_{mach}\right)
\end{equation}
This approximation is independent to the spatial spacing and only depends on the time integration step $\tau$.
This is, to the best of our knowledge, a new state-of-the-art approximation for an incompressible diffeomorphic deformation parameterised by an SVF.
However, it is worth noting that \eqref{eq:jac_error} is only guaranteed
if $m$ remains inside $\mathcal{M}$ during the Euler integration.
Thus, when $\mathcal{M}$ is not equal to the entire spatial domain, this approximation may not hold for large deformations.

\section{Evaluation on synthetic data for incompressible multi-modal registration}

\begin{table}[tb]
	\centering
	\caption{\textbf{Validation on incompressible multi-modal registration.}
	RMSE values between predicted and ground-truth transformations are reported in mm. MAE($\abs{\text{Jacobian} - 1}$) stands for the Mean Absolute Error between the Jacobian map of the predicted transformation and a map uniformly equal to $1$ (perfect incompressibility).
	For both metrics, we reported mean value (standard deviation) over 60 registrations.
	}
	\begin{tabularx}{\textwidth}{c *{4}{Y}}
		\toprule
		 \textbf{Error measures} & Affine & Cubic B-splines & Ours\\
		\midrule
		RMSE (transformation)      & 2,28 (0,68) & 0.96 (0.44) & \textbf{0.90 (0.45)}\\
		MAE($\abs{\text{Jacobian} - 1}$)  & 0.0062 (0.005) & 0.054 (0.02) & \textbf{0.00079 (0.0004)}\\

		\bottomrule
	\end{tabularx}
	\label{tab:multimod_results}
\end{table}

We start by evaluating our method on synthetic incompressible registration in the context of multi-modal MRI data.

\subsubsection{Data generation method} We used T1, T2 and PD brain images from the IXI dataset\footnote{\url{https://brain-development.org/ixi-dataset/}}.
We generated realistic ground-truth incompressible transformations in two steps.
First,
we non-linearly registered a pair of T1 images coming from different patients using a classical cubic B-splines SVF $v_0$.
Second,
we generated a quasi-divergence-free SVF $v$ by projecting
$v_0$
on the space of quasi-divergence-free SVFs.
This corresponds to solving:
\begin{equation}
    \label{eq:data_generation}

        \argmin_{\Theta=\{\phi_i\}_i}
        \frac{1}{2}\norm{v(\Theta) - v_0}^2 \quad
        \text{s.t.}
        \quad \forall i,
        \quad \div{v}(x_i, y_i, z_i) = 0

\end{equation}
where $v$ and $v_0$ are parameterised with classical cubic B-splines as in \eqref{eq:bspline3d}. We note that any potential bias towards the space of divergence-conforming B-splines is avoided in this comparison. We used \texttt{IPOPT}~\cite{Wachter2006} to solve \eqref{eq:data_generation}.

\subsubsection{Evaluation}
We generated a ground-truth quasi-incompressible SVF for 10 patients using the previous generation procedure, taking as fixed images 10 other patients.
Then for a given subject and for any pair of imaging modalities $M_1, M_2 \in \{T1, T2, PD\}$, we warped the first image $I(M_1)$ using the inverse of the ground-truth transformation. The task consists in estimating $v_{GT}$ by registering $I(M_1) \circ \widetilde{\exp}(-v_{GT})$ to $I(M_2)$.
We compared the proposed method with divergence-con\-forming B-splines SVF under incompressibility constraint to a classic registration approach based on cubic B-splines SVF similar to the one used to generate the ground-truth.
We also performed an affine registration to illustrate that the ground-truth transformation is not just affine.

\subsubsection{Implementation details}
We used \texttt{NiftyReg}~\cite{Modat2010} for the diffeomorphic registrations using cubic B-splines SVF.
For both \texttt{NiftyReg} and our incompressible implementation we used the same hyperparameters (grid size 5mm, $3$ levels of pyramid) and objective function, consisting of Normalised Mutual Information (NMI) as a similarity measure (weight $0.95$) and a bending energy regularisation term (weight $0.05$).
The optimiser differs, as we used \texttt{IPOPT} optimiser while a Conjugate Gradient approach is used in \texttt{NiftyReg}.

\subsubsection{Results} \texttt{NiftyReg} and our approach recovered the ground-truth SVF up to a RMSE of the order of the images resolution, as shown in Table.~\ref{tab:multimod_results}. In addition, our method recovered the incompresibility with a higher accuracy.

\begin{figure}[tb!]
    \centering
    \includegraphics[width=\textwidth]{fig1_miccai19_v3_crop.pdf}
    \caption{We aim to register SSFP images at different times of the cardiac cycle.
    During the registration, the myocardium should not compress or expand.
    To impose this constraint we provide a mask of it to our algorithm (dashed in yellow).
    The floating image is warped into $\Phi$(Floating SSFP) so that it matches the reference image, while keeping the deformation of the myocardium incompressible, as shown in the log-Jacobian map.
    Red and blue represent expansion and compression respectively.
    }
    \label{fig:logJacobian}
\end{figure}

\section{Evaluation for myocardial tracking}
We evaluate the proposed method on the STACOM'11 myocardial tracking challenge~\cite{Tobon-Gomez2013}.
In particular, this allows a direct comparison of our framework with iLogDemons\cite{Mansi2011} a state-of-the-art incompressible registration method.

\subsubsection{Data}
The STACOM'11 dataset\footnote{\url{http://stacom.cardiacatlas.org/motion-tracking-challenge/}} contains 4D cine Steady State Free Precession (SSFP) of a full cardiac cycle with $30$ time points for $15$ patients.
Those SSFP come with
$12$ manually tracked landmarks.
We used only 13 of the 15 patients available because we found a shift between the images and the landmarks for two of them which we reported to the organisers.
The coordinates of the landmarks obtained by competitive methods of the challenge are available alongside the original data. This allows a direct comparison with
our method.

\subsubsection{Implementation details}
We used a similar implementation to the one of iLogDemons for this challenge, where we registered the first frame to all subsequent frames for each patient and used a manually delineated mask for the first frame. Yet, we used Local Normalised Cross Correlation as a similarity measure.
We used a bending energy regularisation, and a grid size of $3$mm.

\begin{figure}[bt!]
    \centering
    \begin{subfigure}[t]{0.37\textwidth}
        \centering
        \includegraphics[width=\linewidth]{Figure_boxplot_v2.png}
    \end{subfigure}
    \begin{subfigure}[t]{0.62\textwidth}
        \centering
        \includegraphics[trim={0 0 0 0.8cm},clip, width=\linewidth]{Figure_myocardium_incompressibility_v2.png}
    \end{subfigure}
    \caption{(left) Distances to the landmarks after registering SSFP data.
    Our method achieves similar results to iLogDemons~\cite{Mansi2011} for an incompressible registration.
    UPF~\cite{DeCraene2012} achieves the highest accuracy, but it does not guarantee incompressibility.
    (right) Evolution of the Jacobian values distribution in the incompressible region (myocardium) while registering the first frame to all other time frames. The Jacobian values are uniformly close to 1 for small and moderate deformations of the myocardium. For large deformations (frames at the opposite to the cardiac cycle)
    the mean Jacobian value is close to 1, but the dispersion of the Jacobian value distribution is large.
    This can be attributed to the use of a SVF.
    }
    \label{fig:ssfp_boxplot}
\end{figure}

\subsubsection{Results}
The evaluation of the myocardial tracking is based on the manually annotated landmarks at End Diastole and End Systole.
Using a Wilcoxon signed-rank test, we found that our results, shown in Fig.~\ref{fig:ssfp_boxplot}, are not statistically different from the results of iLogDemons~\cite{Mansi2011} in terms of landmark tracking error.
The log-Jacobian map of Fig.~\ref{fig:logJacobian} illustrates our approximation results \eqref{eq:jac_error}: for a small deformation, our registration framework can guarantee an incompressible deformation in a subregion with high accuracy.
While Fig.~\ref{fig:ssfp_boxplot} illustrates the degradation of this result for larger deformation, i.e. when registering frames at the opposite of the cardiac cycle.

\newpage

\section{Conclusion and Future Work}

\subsubsection{Limitations}
Similar to~\cite{Mang2018,Mansi2011}, the divergence-free constraint is imposed on the stationary velocity field (SVF) in a subregion $\mathcal{M}$, rather than imposing incompressibility on $\mathcal{M}$.
A voxel that is originally in $\mathcal{M}$ might exit $\mathcal{M}$ during the integration of the SVF and start being transported by velocity vectors that are not divergence-free (see supplemental material for a mathematical justification).
In this case, the deformation of this voxel is no longer incompressible.
As a result, using an SVF for incompressible registration applies only for deformations that are small enough or when the whole spatial domain is constrained to be divergence-free.
Investigating other diffeomorphic parameterisations is left for future work.

\subsubsection{Advantages compared to previous methods}
Our method for incompressible registration relies on the parametrisation of the velocity field by a divergence conforming B-splines~\cite{Evans2013a}.
In contrast to classical B-splines used in~\cite{Mang2018,Mansi2011}, it guarantees that the divergence of the velocity field is still a B-spline.
This parametrisation along with constrained optimisation methods allows us to impose the velocity field to be divergence-free up to machine precision ($10^{-16}$ in our experiments). This is irrespective of the grid resolution chosen for the velocity field.
As a result, the proposed method
is scalable to 3D images with high resolution.

We also proved an error bound for the incompressibility of the deformation for the proposed method in the case an Euler integration of the velocity field is used \eqref{eq:jac_error}.
The study of the error bound for other integrators, that may require additional interpolations (e.g. scaling-and-squaring), is left for future works.
Additionally, previous (quasi-)incompressible registration methods~\cite{Mang2018,Mansi2011} have been limited to using SSD as image similarity metric. We are proposing the first incompressible non-linear registration method that supports any smooth image similarity measure and spatial regularization.

\subsubsection*{Acknowledgments}
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk{\l}odowska-Curie grant agreement TRABIT No 765148;
Wellcome [203148/Z/16/Z; 203145Z/16/Z; WT101957], EPSRC [NS/A000049/1; NS/A000050/1; NS/A000027/1; EP/L016478/1].
TV is supported by a Medtronic / RAEng Research Chair [RCSRF1819\textbackslash7\textbackslash34].

\bibliographystyle{splncs04.bst}