Small sphere distributions for directional data with application to medical imaging

TL;DR

This paper introduces new small-sphere distribution models for multivariate directional data, particularly useful in medical imaging for analyzing shape and shape changes of 3D objects with directions concentrated on small circles.

Contribution

The paper develops novel multivariate small-sphere distributions that model association among directions and facilitate inference, with applications in medical imaging and shape analysis.

Findings

Effective modeling of shape changes in 3D objects.

Application to human knee gait data.

Advantages over existing small-circle data analysis methods.

Abstract

We propose new small-sphere distributional families for modeling multivariate directional data on $(\mathbb{S}^{p-1})^K$ for $p \ge 3$ and $K \ge 1$. In a special case of univariate directions in $\Re^3$, the new densities model random directions on $\mathbb{S}^2$ with a tendency to vary along a small circle on the sphere, and with a unique mode on the small circle. The proposed multivariate densities enable us to model association among multivariate directions, and are useful in medical imaging, where multivariate directions are used to represent shape and shape changes of 3-dimensional objects. When the underlying objects are rotationally deformed under noise, for instance, twisted and/or bend, corresponding directions tend to follow the proposed small-sphere distributions. The proposed models have several advantages over other methods analyzing small-circle-concentrated data, including inference procedures on the association and small-circle fitting. We demonstrate the use of the proposed multivariate small-sphere distributions in analyses of skeletally-represented object shapes and human knee gait data.

Tables

Table 1. Table 1 : Small-sphere distributions (top three rows) and methods (bottom three rows) developed in this paper. Items with “ × \times ” mark are beyond the scope of this paper.

	First kind (S1)	Second kind (S2)
Univariate	S1	S2
Multivariate (indep.)	iMS1	iMS2
Multivariate (dep.)	$\times$	GMS2	MS2
Simulation	Gibbs sampling	$\times$	Exact sampling
Estimation	Approximate m.l.e.	$\times$	Approximate m.l.e.
Hypothesis testing	Likelihood ratio	$\times$	Likelihood ratio

Table 2. Table 2 : Small-circle estimation performances for univariate data on 𝕊 2 superscript 𝕊 2 \mathbb{S}^{2} from Fig.  4 . Means (standard deviations) of the angular product errors in degrees ( 32 ) are shown.

Table 3. Table 3 : Small-circles estimation performances for bivariate data on 𝕊 2 superscript 𝕊 2 \mathbb{S}^{2} from Fig. A6 in the online supplementary material. Means (standard deviations) of the angular product errors in degrees ( 32 ) are shown.

Table 4. Table 4 : Concentration and association parameter estimates for bivariate data on ( 𝕊 2 ) 2 superscript superscript 𝕊 2 2 (\mathbb{S}^{2})^{2} from Fig. A6 in the online supplmentary material. Means (standard deviations) of the estimates (from 100 repetitions). The column headings show the true parameters.

Equations

$$\mathcal{C}(\mu ,\nu )=\{x\in {\mathbb{S}}^{p-1}|\delta (\mu ,x)=\arccos (\nu )\},\nu \in (-1,1).$$

$$f_{\text{vMF}}(x;\mu ,\kappa )={\left(\frac{\kappa }{2}\right)}^{p/2-1}\frac{1}{\Gamma (p/2)I_{p/2-1}(\kappa )}\exp \{\kappa {\mu }^{\top }x\},x\in {\mathbb{S}}^{p-1}.$$

$$f_{\text{BM}}(x;\mu ,\kappa ,\nu )=\frac{1}{\alpha (\kappa ,\nu )}\exp \{-\kappa ({\mu }^{\top }x-\nu {)}^2\},x\in {\mathbb{S}}^{p-1},$$

$$g_{\text{vMF}}(s,\varphi ;\kappa )$$

$$g_{\text{BM}}(s,\varphi ;\kappa ,\nu )$$

$$f_{\text{S1}}(x;{\mu }_0,{\mu }_1,{\kappa }_0,{\kappa }_1)$$

$$f_{\text{S2}}(x;{\mu }_0,{\mu }_1,{\kappa }_0,{\kappa }_1)$$

$$y:=(x_2,\dots ,x_p)/\sqrt{1-s^2}\in {\mathbb{S}}^{p-2},$$

$$g_{\text{S1}}(s,y;{\mu }_1,{\kappa }_0,{\kappa }_1)$$

$$g_{\text{S2}}(s,y;{\mu }_1,{\kappa }_0,{\kappa }_1)$$

$$f_{S1}(x;\gamma ,A)=\frac{1}{\alpha (\gamma ,A)}\exp \{{\gamma }^{\top }x-x^{\top }Ax\},$$

$$\alpha (\gamma ,A)=a({\kappa }_0,{\kappa }_1,\nu )\exp \left\{{\kappa }_0{\nu }^2\right\}.$$

$$f_{iMS1}(\text{x})\propto \exp \left\{{\Gamma }^{\top }\text{x}-{\text{x}}^{\top }\text{A}\text{x}\right\}.$$

$$g_{iMS2}(\text{s},\text{y})\propto \exp \left\{H^{\top }\text{s}-{\text{s}}^{\top }{D}_{{\kappa }_0}\text{s}+{\text{M}}^{\top }\mathrm{vec}(\text{y})\right\},$$

$$g_{\text{GMS2}}(\text{s},\text{y};H,D_{{\kappa }_0},\text{M},\text{B})$$

$$B^{k,l}$$

$$\widetilde{{\mu }_{2(k)}}=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\widetilde{{\mu }_{1(k)}}.$$

$${\widetilde{{\mu }_{1(k)}}}^{\top }{y}_{(k)}=\cos ({\varphi }_{(k)}-{\zeta }_{(k)}),{\widetilde{{\mu }_{2(k)}}}^{\top }{y}_{(k)}=\sin ({\varphi }_{(k)}-{\zeta }_{(k)}).$$

$$c(\mathbf{\varphi },\mathbf{\zeta })$$

$$s(\mathbf{\varphi },\mathbf{\zeta })$$

$$g_{\text{MS2}}(\text{s},\mathbf{\varphi };H,D_{{\kappa }_0},{\mathbf{\kappa }}_1,\mathbf{\zeta },\mathbf{\Lambda })$$

$${\mathbf{\kappa }}_1^{\top }c(\mathbf{\varphi },\mathbf{\zeta })+\frac{1}{2}s(\mathbf{\varphi },\mathbf{\zeta }{)}^{\top }\mathbf{\Lambda }s(\mathbf{\varphi },\mathbf{\zeta })$$

$$a(\kappa_{0},\kappa_{1},\nu)=2\pi^{p/2}|\mbox{\boldmath{$\Psi$}}|^{-1/2}g(1)\exp\left(\xi^{T}\Psi\xi+h-\kappa_{0}\nu^{2}\right).$$

$$K_g^{(j)}(t)=\frac{(j-1)!}{2}\left(\frac{1}{({\kappa }_0+h-t{)}^j}+\frac{p-1}{(h-t{)}^j}\right)+\frac{j!}{4}\left(\frac{{\nu }^2(2{\kappa }_0+{\kappa }_1{)}^2}{({\kappa }_0+h-t{)}^{j+1}}+\frac{{\kappa }_1^2(1-{\nu }^2)}{(h-t{)}^{j+1}}\right).$$

$$\hat{g}(1)=(2\pi K_g^{(2)}(\hat{t}){)}^{-1/2}\exp (K_g(\hat{t})-\hat{t}+T),$$

$$\underset{{\mu }_0\in {\mathbb{S}}^{p-1}}{\min }\left[\frac{1}{n}\Vert {\mathbb{X}}^{\top }{\mu }_0-\nu 1_n{\Vert }^2-\lambda ({\mu }_0^{\top }{\mu }_0-1)\right],$$

$${\hat{\mu }}_0=\nu (S-\hat{\lambda }{I}_p{)}^{-1}\bar{x}.$$

$${\mu }_1=\cos (\phi ){\mu }_0+\sin (\phi ){\gamma }^{\ast }.$$

$${\tilde{\ell }}_{{\mu }_0}(\phi ,{\kappa }_0,{\kappa }_1)=-n\log \hat{a}({\kappa }_0,{\kappa }_1,\phi )-{\kappa }_0\sum _{i=1}^n({\mu }_0^{\top }{x}_i-\cos \phi {)}^2+{\kappa }_1\sum _{i=1}^n{x}_i^{\top }(\cos (\phi ){\mu }_0+\sin (\phi ){\hat{\gamma }}^{\ast }).$$

$$\underset{{\mu }_0\in {\mathbb{S}}^{p-1}}{\min }\left[\frac{1}{n}\sum _{j=1}^K({\kappa }_{0(j)}\Vert {\mathbb{X}}_{(j)}^{\top }{\mu }_0-{\nu }_{(j)}{1}_n{\Vert }^2)-\lambda ({\mu }_0^{\top }{\mu }_0-1)\right],$$

$${\hat{\mu }}_0={\mathrm{argmin}}_{{\mu }_0}{\ell }_n({\mu }_0,{\hat{\eta }}_{{\mu }_0}),$$

$${\hat{\eta }}_{{\mu }_0}={\mathrm{argmin}}_{\eta }{\ell }_n({\mu }_0,\eta )$$

Full text

Keywords: Bingham-Mardia distribution, directional data, skeletal representation, small circle, small sphere, likelihood ratio test, maximum likelihood estimation, von Mises-Fisher distribution.

1 Introduction

In medical imaging, accurately assessing and correctly diagnosing shape changes of internal organs is a major objective of a substantial challenge. Shape deformations can occur through long-term growth or necrosis as well as by short-term natural deformations. In view of surgery and radiation therapy, it is important to model all possible variations of object deformations by both long- and short-term changes, in order to control the object’s exact status and shape at treatment time. Rotational deformations such as rotation, bending, and twisting form a key sub-category of possible shape changes. For instance, shape changes of hippocampi in the human brain have been shown to mainly occur in the way of bending and twisting (Joshi et al., 2002; Pizer et al., 2013).

For the task of modeling 3D objects an abundance of approaches have been introduced. Closely related to our work are landmark-based shape models (Cootes et al., 1992; Dryden and Mardia, 1998; Kurtek et al., 2011) where a solid object is modeled by the positions of surface points, chosen either anatomically, mathematically or randomly. A richer family of models is obtained by attaching directions normal to the sampled surface points. More generally, in skeletal representations (called s-reps, Siddiqi and Pizer, 2008), an object is modeled by the combination of skeletal positions (lying on a medial sheet inside of the object) and spoke vectors (connecting the skeletal positions with the boundary of the object). In these models, describing the variation of rotational deformations can be transformed into a problem of exploring the motion of directional vectors on the unit two-sphere. As argued in Schulz et al. (2015), directional vectors representing rotational deformations tend to be concentrated on small circles on the unit sphere; a toy data example in Fig. 1 shows a typical pattern of such observations.

Motivated by the analysis of such s-rep data, we propose new distributional families and their multivariate extensions in order to model such directional data on the unit sphere $\mathbb{S}^{p-1}=\{x\in\Re^{p}\,|\,\|x\|=1\}$ in arbitrary dimension $p\geq 3$. (Throughout the paper, $\|x\|=(x^{\top}x)^{1/2}$ is the usual 2-norm of the vector $x$.) To precisely describe the targeted data situation, we define a $(p-2)$-dimensional subsphere of $\mathbb{S}^{p-1}$ as the set of all points equidistant from $\mu\in\mathbb{S}^{p-1}$, denoted by

[TABLE]

Here, $\delta(u,v)=\arccos(u^{\top}v)$ is the geodesic distance between $u,v\in\mathbb{S}^{p-1}$. The subsphere is called a great subsphere if $\nu=0$ and a proper small subsphere if $\nu\neq 0$. Note that $\mathcal{C}(\mu,\nu)\subset\mathbb{S}^{p-1}$ is well-defined for all $p>1$. For the special case of $p=3$, $\mathcal{C}(\mu,\nu)$ is a circle, a great circle if $\nu=0$, and a proper small circle if $\nu\neq 0$. To model the data in Fig. 1, one may naively use the Bingham-Mardia distribution (called “BM” hereafter, Bingham and Mardia, 1978), which is a family of densities on $\mathbb{S}^{2}$ with a modal ridge along a small circle. However, typical observations we encountered in applications do not uniformly spread over the full circle, and the BM distribution does not fit well, as shown in Fig. 1(b). Moreover, when by a single observation multiple directional vectors are provided, that is, data are on a polysphere $(\mathbb{S}^{2})^{K}$, to the knowledge of the authors, there is no tool available to date, to model dependencies between directions.

In this paper, we propose two types of new distributional families for random directional vectors on $\mathbb{S}^{p-1}$, which we call small-sphere distributions of the first (S1) and second (S2) kind. If $p=3$, the proposed distributions may be called small-circle distributions. These two distributional families are designed to have higher densities on $\mathcal{C}(\mu,\nu)$ and to have a unique mode on $\mathcal{C}(\mu,\nu)$. An example of a small-sphere density, fitted to the toy data is shown in Fig. 1(c). The new densities are natural extensions of the BM distribution with an additional term explaining a decay from a mode. If the additional term is a von Mises-Fisher (vMF) density on $\mathbb{S}^{p-1}$, we obtain the S1, which is a subfamily of the general Fisher-Bingham distribution (Mardia, 1975; Mardia and Jupp, 2000). On the other hand, if the additional term is a vMF density on the small sphere ($\cong\mathbb{S}^{p-2}$), we obtain the S2 distribution, in which case the “horizontal” and “vertical” components of the directional vectors are independent of each other.

Several multivariate extensions of the new distributions to $(\mathbb{S}^{p-1})^{K}$, $K\geq 2$, are discussed as well. In particular, we show that a special case, called MS2, of our multivariate extensions is capable of modeling dependent random vectors. It has a straightforward interpretation, and we provide for fast estimation of its parameters. This MS2 distribution is specifically designed with s-rep applications in mind. In particular, s-rep data from rotationally-deformed objects have directional vectors that are “rotated together,” share a common axis of rotation, and are “horizontally dependent” (when the axis is considered to be vertically positioned). The component-wise independence of the S2 distributions plays a key role in this simple and interpretable extension. We discuss here likelihood-based parameter estimation and testing procedures of the multivariate distributions.

The contribution of this paper is summarized in Table 1. While the new distributions contribute to the literature of directional distributions (Mardia and Jupp, 2000), the proposed estimation procedures for the S1, S2 and MS2 parameters can be thought of as a method of fitting small-subspheres to data, which has been of separate interest. Nonparametric least-squares type solutions for such problem dates back to Mardia and Gadsden (1977), Gray et al. (1980), and Rivest (1999). Jung et al. (2012) proposed to recursively fitting small-subspheres in dimension reduction of directional and shape data. Pizer et al. (2013) proposed to combine separate small-circle fitting results in the analysis of s-rep data. In a similar spirit, Jung et al. (2011) and Schulz et al. (2015) also considered fitting small-circles in applications to s-rep analysis. In a simulation study, we show that our estimators provide smaller mean angular errors in small-circle fits than recent developments listed above. Moreover, our parametric framework provides a formal large-sample likelihood ratio test procedure, which is applicable to a number of hypothesis settings and is more powerful in testing dependence among multivariate directions than competing methods.

The rest of the article is organized as follows. In Section 2, we introduce the proposed densities of the S1 and S2 distributions and discuss their multivariate extensions including the MS2 distribution. Procedures of obtaining random variates from the proposed distributions are also proposed and discussed. In Section 3, algorithms to obtain approximate maximum likelihood estimators of the parameters are proposed and discussed. In Section 4, we introduce several hypotheses of interest and procedures of large-sample approximate likelihood-ratio tests. In Section 5, we empirically show that the proposed models are superior over other methods for small-circle estimation and for estimating dependencies in the multivariate setting, and show that the proposed testing procedures effectively prevent overfitting. In Sections 6 and 7 we demonstrate applications of the new multivariate distributions to analyze models that represent human organs and knee motions. The appendix and the online supplementary material contain technical details and additional numerical results.

2 Parametric small-sphere models

First we introduce two classical spherical densities, then we suitably combine them for our purposes.

2.1 Two classical distributions on $\mathbb{S}^{p-1}$

The von Mises-Fisher (vMF) distribution (Mardia and Jupp, 2000, p.168) is a fundamental unimodal and isotropic distribution for directions with density

[TABLE]

Here, $\Gamma$ is the gamma function and $I_{v}$ is the modified Bessel function of the first kind and order $v$. The parameter $\mu\in\mathbb{S}^{p-1}$ locates the unique mode with $\kappa\geq 0$ representing the degree of concentration.

The Bingham-Mardia (BM) distribution was introduced by Bingham and Mardia (1978) to fit data in $\mathbb{S}^{2}$ that cluster near a small circle $\mathcal{C}(\mu,\nu)$. For an arbitrary dimension $p\geq 3$, the BM density is given by

[TABLE]

where $\alpha(\kappa,\nu)>0$ is the normalizing constant.

For our purpose of generalizing these distributions, we represent the variable $x\in\mathbb{S}^{p-1}$, $p\geq 3$, by spherical angles $\phi_{1},\ldots,\phi_{p-1}$ satisfying $\cos\phi_{1}=\mu^{\top}x$. Setting $s:=\cos\phi_{1}\in[-1,1]$ and $\phi:=(\phi_{2},\ldots,\phi_{p-1})\in[0,\pi]^{p-3}\times[0,2\pi)$, the random vector $(s,\phi)$ following the von Mises-Fisher (1) or Bingham-Mardia (2) distribution has the respective density:

[TABLE]

In consequence, for both distributions, $s$ and $\phi$ are independent, and the marginal distribution of $\phi$, which parametrizes a co-dimension 1 unit sphere $\mathbb{S}^{p-2}$, is uniform. In (3), the marginal distribution of $s$ is a shifted exponential distribution truncated to $s\in[-1,1]$, while in (4) the marginal distribution of $s$ is a normal distribution truncated to $s\in[-1,1]$. Both densities are isotropic, i.e. rotationally symmetric with respect to $\mu$. The vMF density is maximal at the mode $\mu$ and decreases as the latitude $\phi_{1}$ increases, while the BM density is uniformly maximal on the small-sphere $\mathcal{C}(\mu,\nu)$ and decreases as $\phi_{1}$ deviates from $\arccos(\nu)$.

2.2 Small-sphere distributions of the first and second kind

The proposed small-sphere densities of the first and second kind on $\mathbb{S}^{p-1}$, for $x=(x_{1},\ldots,x_{p})\in\mathbb{S}^{p-1}$ with parameters $\mu_{0},\mu_{1}\in\mathbb{S}^{p-1}$, $\nu=\mu_{0}^{\top}\mu_{1}\in(-1,1)$, $\kappa_{0}>0$, $\kappa_{1}>0$, are given by

[TABLE]

respectively, where $a(\kappa_{0},\kappa_{1},\nu)$ and $b(\kappa_{0},\kappa_{1},\nu)$ are normalizing constants. Here, $P_{\mu_{0}}$ denotes the matrix of orthogonal projection to the orthogonal complement of $\mu_{0}$; $P_{\mu_{0}}=I_{p}-\mu_{0}\mu_{0}^{\top}$, where $I_{p}$ is the identity matrix. (In (6), we use the convention $0/0=0$.)

These distributions are well-suited to model observations that are concentrated near the small sphere $\mathcal{C}(\mu_{0},\nu)$ but are not rotationally symmetric. The first kind (5) is a natural combination of the vMF (1) and BM (2) distributions. The parameter $\mu_{0}$ represents the axis of the small sphere $\mathcal{C}(\mu_{0},\nu)$, while $\mu_{1}$ gives the mode of the distribution, which, by the definition of $\nu$, is on the small sphere $\mathcal{C}(\mu_{0},\nu)$. These parameters, $\mu_{0},\mu_{1}$, $\nu$, are illustrated in Fig. 1(a) for the $p=3$ case. The parameter $\kappa_{0}$ controls the vertical concentration towards the small sphere (with an understanding that $\mu_{0}$ is arranged vertically). In (5), $\kappa_{1}$ controls the isotropic part of the concentration around the mode, forcing the density to decay from $\mu_{1}$.

The rationale for the second kind (6) is better understood with a change of variables. Let us assume for now that $\mu_{0}=(1,0,\ldots,0)^{\top}$. For any $x=(x_{1},\ldots,x_{p})^{\top}\in\mathbb{S}^{p-1}$, write $s:=x_{1}=\mu_{0}^{\top}x$. If the spherical coordinate system $(\phi_{1},\ldots,\phi_{p-1})$ as defined for (4) is used, then $s=\cos\phi_{1}$. The “orthogonal complement” of $s$ is denoted by

[TABLE]

where the vector $y$ is obtained from the relation $P_{\mu_{0}}x/\|P_{\mu_{0}}x\|=(0,y)\in\mathbb{S}^{p-1}$. Similarly, define $\widetilde{\mu_{1}}\in\mathbb{S}^{p-2}$ as the last $p-1$ coordinates of $P_{\mu_{0}}\mu_{1}/\|P_{\mu_{0}}\mu_{1}\|$. Then the random vector $(s,y)\in[-1,1]\times\mathbb{S}^{p-2}$ from the S1 or S2 has densities

[TABLE]

respectively, for $s\in[-1,1],\,y\in\mathbb{S}^{p-2}$. The subtle difference is that for the first kind (8), the “vMF part” (the second term in the exponent) is not statistically independent from the “BM part”, while it is true for the second kind (9). That is, $s$ and $y$ are independent only in the second kind. Accordingly, in (9), $\kappa_{1}$ controls the horizontal concentration towards the mode $\mu_{1}$. The parameters $\mu_{0},\mu_{1}$ and $\kappa_{0}$ of the second kind have the same interpretations as those of the first kind.

We use the notation $X\sim\textrm{S1}(\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$ and $Y\sim\textrm{S2}(\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$ for random directions $X,Y\in\mathbb{S}^{p-1}$ following small-sphere distributions of the first and second kind with parameters $(\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$, respectively. The proposed distributions are quite flexible and can fit a wide range of data. In Figure 2, we illustrate the S1 densities (5) with various values of the concentration parameters $\kappa_{0},\kappa_{1}$. In all cases, the density is relatively high near the small circle $\mathcal{C}(\mu_{0},\nu)$ and has the mode at $\mu_{1}\in\mathcal{C}(\mu,\nu)$. Despite the difference in their formulations, the S2 densities (6) look similar to S1 densities for each fixed parameter-set. We refer to the online supplementary material for several visual examples of the S2 density.

The following invariance properties are proven in the Appendix.

Proposition 1**.**

Let $X,Y\in\mathbb{S}^{p-1}$ be random directions with $X\sim\textrm{S1}(\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$ and $Y\sim\textrm{S2}(\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$ and let $B$ be a full rank $p\times p$ matrix.

• (i)

If $B\mu_{0}=\mu_{0}$ and $B\mu_{1}=\mu_{1}$ then $X$ and $BX$ have the same distribution and so do $Y$ and $BY$.

• (ii)

If $B$ is orthogonal and $X,BX$ have the same distribution (or $Y,BY$ have the same distribution), then $B\mu_{0}=\mu_{0}$ and $B\mu_{1}=\mu_{1}$.

• (iii)

$X\sim\textrm{S1}(-\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$* and $Y\sim\textrm{S2}(-\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$.*

An example for the matrix $B$ in Proposition 1(i) is the reflection matrix $B=I_{p}-2UU^{\top}$, where $U=[u_{3},\ldots,u_{p}]$ is such that $[u_{1},\ldots,u_{p}]$ is a $p\times p$ orthogonal matrix whose column vectors $u_{1}$ and $u_{2}$ generate $\mu_{0}$ and $\mu_{1}$.

Remark 1*.*

The S1 distribution is a special case of the Fisher-Bingham distribution (Mardia, 1975). Following the notation of Kent (1982), the S1 distribution may be labeled as a FB6 distribution, in the special case of $p=3$, emphasizing the 6-dimensional parameter space. In terms of the general parameterization of the Fisher-Bingham density (cf. Mardia and Jupp, 2000, p.174), we write $\gamma=2\kappa_{0}\nu\mu_{0}+\kappa_{1}\mu_{1}$ and $A=\kappa_{0}\mu_{0}\mu_{0}^{\top}$, so that the S1 density (5) is expressed as

[TABLE]

where

[TABLE]

This relation to the general Fisher-Bingham distribution facilitates random data generation and maximum likelihood estimation, shown later in Sections 2.4 and 3.1.

2.3 Multivariate extensions

The univariate small-sphere distributions (5–6) are now extended to model a tuple of associated random directions, $\textbf{X}=(X_{(1)},\ldots,X_{(K)})\in(\mathbb{S}^{p-1})^{K}$. We confine ourselves to a special case where the marginal distributions of $X_{(k)}$ have a common “axis” parameter $\mu_{0}$, but relaxing this condition is straightforward. We begin by introducing multivariate small-sphere distributions for independent random directions, denoted by iMS1 and iMS2.

Independent extensions.

Suppose that, in the $K$-tuple of random directions X, each $X_{(k)}\in\mathbb{S}^{p-1}$ is marginally distributed as S1$(\mu_{0},\mu_{1(k)},\kappa_{0(k)},\kappa_{1(k)})$. Throughout, we assume that $\nu_{(k)}=\mu_{0}^{\top}\mu_{1(k)}\in(-1,1)$ so that the underlying small spheres do not degenerate. If the components of X are mutually independent, then the joint density evaluated at $\textbf{x}\in(\mathbb{S}^{p-1})^{K}$ is

[TABLE]

Here, $\Gamma=[\gamma_{(1)},\ldots,\gamma_{(K)}]^{\top}$, where $\gamma_{(k)}=2\kappa_{0(k)}\nu_{(k)}\mu_{0}+\kappa_{1(k)}\mu_{1(k)}$, and $\textbf{A}=D_{\kappa_{0}}\otimes(\mu_{0}\mu_{0}^{\top})$, where $D_{\kappa_{0}}=\textrm{diag}(\kappa_{0(1)},\ldots,\kappa_{0(K)})$. Each marginal density is of the form (10).

If each component is marginally distributed as S2$(\mu_{0},\mu_{1(k)},\kappa_{0(k)},\kappa_{1(k)})$, then writing the density in terms of $(s,y)$ as done for (9) facilitates our discussion. First, we decompose each $x_{(k)}$ into $s_{(k)}=\mu_{0}^{\top}x_{(k)}\in[-1,1]$ and $y_{(k)}\in\mathbb{S}^{p-2}$ as defined in (7). Further, we denote by $\widetilde{\mu_{1(k)}}$ the scaled projection of $\mu_{1(k)}$ as done for the univariate case. Then an independent multivariate extension for the S2 model can be expressed as the joint density of $\textbf{s}=(s_{(1)},\ldots,s_{(K)})$ and $\textbf{y}=(y_{(1)},\ldots,y_{(K)})$,

[TABLE]

where $H=(2\kappa_{0(1)}\nu_{(1)},\ldots,2\kappa_{0(K)}\nu_{(K)})$ and $\textbf{M}=\mathrm{vec}(\kappa_{1(1)}\widetilde{\mu_{1(1)}},\ldots,\kappa_{1(K)}\widetilde{\mu_{1(K)}})$ while $\mathrm{vec}(\cdot)$ denotes the column-wise vectorization of a matrix.

Vertical and horizontal dependence.

Based on (12) and (13), we now contemplate on dependent models. Obviously, if we allow in (12) nonzero offdiagonal entries of A, then we obtain a dependent modification of the S1 model. With our applications in mind, however, we aim at modeling a specific structure of dependence that is natural to the variables $(\textbf{s},\textbf{y})$ in (13).

If $s_{(1)},\ldots,s_{(K)}$ are dependent, we speak of vertical dependence; if $y_{(1)},\ldots,y_{(K)}$ are dependent, we speak of horizontal dependence. In practice, when we deal with small-circle concentrated directional data, association among these vectors usually occurs along small-circles with independent vertical errors. For example, when a 3D object is modeled by skeletal representations, as described in more detail in Section 6 and visualized in Fig. 6, a deformation of the object is measured by the movements of directional vectors on $\mathbb{S}^{2}$. When a single rotational deformation (such as bending, twisting or rotation) occurs, all the directions move along small-circles with a common axis $\mu_{0}$. In this situation, the longitudinal variations along the circles are dependent of each other because nearby spoke vectors are under the effect of similar deformations. (Examples of such longitudinal dependencies can be found in Section 6 as well as in Schulz et al. (2015).) Adding such a horizontal (or longitudinal) dependence to a multivariate S1 model requires a careful introduction and parametrization of the offdiagonal entries of A in (12). This is not straightforward, and we leave it for future work. On the other hand, it is feasible to extend the S2 model by generalizing the “vMF part” of y, the last term in the exponent of (13), to a Fisher-Bingham type.

To this end, we introduce a parameter matrix B to model general quadratics in $\mathrm{vec}(\textbf{y})$. This allows to write the densities for a general multivariate small-sphere distribution of the second kind (GMS2) as follows:

[TABLE]

where $H,D_{\kappa_{0}}$ and M as defined in (13), and $T_{1}(H,D_{\kappa_{0}})$ and $T_{2}(\textbf{M},\textbf{B})$ are normalizing constants. We set $\textbf{B}=(B_{k,l})_{k,l=1}^{K}$, $B_{k,l}=(b^{(k,l)}_{i,j})_{i,j=1}^{p-1},$ as a block matrix with vanishing blocks $B_{k,k}=0$ on the diagonal. The submatrix $B_{k,l}$ models the horizontal association between $y_{(k)}$ and $y_{(l)}$. The fact that $z^{\top}\textbf{B}z=z^{\top}\textbf{B}^{\top}z=\frac{1}{2}z^{\top}(\textbf{B}+\textbf{B}^{\top})z$ for any vector $z\in\Re^{(p-1)K}$ allows us to assume without loss of generality that B is symmetric.

An MS2 distribution on $(\mathbb{S}^{2})^{K}$.

As a viable submodel for the practically important case $p=3$, we propose to use a special form for the offdiagonal blocks $B_{k,l}$ of B. In particular, with $\lambda_{k,l}$ representing the degrees of association between $y_{(k)}$ and $y_{(l)}$, we set

[TABLE]

where $\begin{pmatrix}\widetilde{\mu_{1(k)}}&\widetilde{\mu_{2(k)}}\\ \end{pmatrix}$ is the rotation matrix given by setting

[TABLE]

The density (14) with the above parameterization of B will be referred to as a multivariate S2 distribution (MS2) for data on $(\mathbb{S}^{2})^{K}$; its angular representation will be derived in (16) below.

Our choice of the simple parametrization (2.3) does not restrict the modeling capability of the general model (14), and has some advantages in parameter interpretations and also in estimation. To see this, we resort to use an angular representation for y (available to this $p=3$ case). For each $k$, define $\phi_{(k)}$ and $\zeta_{(k)}$ such that $y_{(k)}=(\cos\phi_{(k)},\sin\phi_{(k)})^{\top}$ and $\widetilde{\mu_{1(k)}}=(\cos\zeta_{(k)},\sin\zeta_{(k)})^{\top}$. Accordingly, the inner products appearing in (14) can be expressed as

[TABLE]

Let $\boldsymbol{\phi}=(\phi_{(1)},\ldots,\phi_{(K)})^{\top}$, $\boldsymbol{\zeta}=(\zeta_{(1)},\ldots,\zeta_{(K)})^{\top}$, $\boldsymbol{\kappa}_{1}=(\kappa_{1(1)},\ldots,\kappa_{1(K)})^{\top}$,

[TABLE]

and $\boldsymbol{\Lambda}=(\lambda_{k,l})_{k,l=1}^{K}$ where $\lambda_{k,l}(=\lambda_{l,k})$ for $k\neq l$ is the association parameter used in (2.3), and $\lambda_{k,k}$ is set to zero. The density of the MS2 distribution, in terms of $(\textbf{s},\boldsymbol{\phi})$, is then

[TABLE]

From (16), it can be easily seen that the “horizontal angles” $\boldsymbol{\phi}$ follow the multivariate von Mises distribution (Mardia et al., 2008) and are independent of the vertical component s. As we will see later in Section 3.2, this facilitates estimation for the MS2 distributions. Moreover, since

[TABLE]

for large enough concentrations, $\boldsymbol{\phi}$ is approximately multivariate normal with mean $\boldsymbol{\zeta}$ and precision matrix $\Sigma^{-1}$, where $(\Sigma^{-1})_{kk}=\kappa_{1(k)}$ and $(\Sigma^{-1})_{kl}=-\lambda_{kl}$ for $1\leq k\neq l\leq K$. These parameters are naturally interpreted as partial variances and correlations. This interpretation of the parameters as entries of a precision matrix is most immediate under the MS2, but is not under the general case.

2.4 Random data generation

Generating pseudo-random samples from the S1 and S2 distributions are important in simulations and in developments of computer-intensive inference procedures.

For simulation of the S1 (5) and iMS1 (12) distribution, the fact that each marginal distribution of the iMS1 is a special case of the Fisher-Bingham is handy. Thereby, one can use the Gibbs sampling procedure developed for generating Fisher-Bingham-variate samples Hoff (2009).

For simulation of the S2 (6), iMS2 (13), and MS2 (16) distribution, we take advantage of the independence between the pair ($\textbf{s},\textbf{y}$). As we assume vertical independence (i.e., $s_{(1)},\ldots,s_{(K)}$ are independent), each $s_{(k)}$ can be sampled separately. Therefore, sampling from the MS2 distribution amounts to independently drawing samples from a truncated normal distribution (for $s_{(k)}$) and from a multivariate von Mises distribution (for y). Specifically, to sample $\mathbf{x}=(x_{(1)},\ldots,x_{K)})$ from $\textrm{MS}2(\mu_{0},\boldsymbol{\mu}_{1},\boldsymbol{\kappa}_{0},\boldsymbol{\kappa}_{1},\boldsymbol{\Lambda})$, the following procedure can be used.

Step 1. For each $k$, sample $s_{(k)}$ from the truncated normal distribution with mean $\nu_{(k)}$ and variance $1/(2\kappa_{0(k)})$, truncated to the interval (-1,1).

Step 2. For the S2 or iMS2 model, sample each $y_{(k)}\in\mathbb{S}^{p-2}$ in $\textbf{y}=(y_{(1)},\ldots,y_{(K)})$ independently from the von Mises distribution with mean $(1,0,\ldots,0)$ and concentration $\kappa_{1(k)}$; for the MS2 distribution (when $p=3$), sample the $K$-tuple $\textbf{y}\in(\mathbb{S}^{1})^{K}$ directly from the multivariate von-Mises distribution with mean $(1,0)$ and precision parameters $\boldsymbol{\kappa}_{1}$ and $\boldsymbol{\Lambda}$.

Step 3. For each $k$, let $E_{(k)}$ be a $p\times p$ orthogonal matrix with $(\mu_{0},P_{\mu_{0}}\mu_{1(k)}/\|P_{\mu_{0}}\mu_{1(k)}\|)$ being the first two column vectors. Set $x_{(k)}=E_{(k)}^{\top}\left(s_{(k)},(1-s_{(k)}^{2})^{-1/2}y_{(k)}\right)$.

In our experiments, sampling from the S2 and MS2 distributions is much faster than from the S1. In particular, when the dimension $p$ or the concentration level is high, the Markov chain simulations for the S1 appear to be sluggish. Some examples of random samples from the S2, iMS2 and MS2 distributions are shown in Fig. 3. The small-circles $C(\mu_{0},\nu_{(k)})$ are also overlaid in the figure. Notably, the MS2 sample in the rightmost panel clearly show a horizontal dependence.

3 Estimation

3.1 Approximate maximum likelihood estimation for S1 and iMS1

The standard way to estimate parameters of the S1 is to use the maximum likelihood estimates (m.l.e). However, it does not seem possible to obtain explicit expressions of the m.l.e., partly due to having no closed-form expression of the normalizing constant (5). We propose to use an approximate m.l.e., obtained by iterating between updates for $\mu_{0}$ and for $(\mu_{1},\kappa_{0},\kappa_{1})$, where each separate problem can be solved efficiently. The proposed estimation procedure, given below, only needs to specify an initial value for $\nu=\mu_{0}^{\top}\mu_{1}\in(-1,1)$. Our procedure naturally extends to estimation for the iMS1 distribution, which will also be discussed.

As a preparation, we first describe an approximation of the normalizing constant. For this, we follow Kume and Wood (2005), who used saddle-point density approximations for approximating normalizing constants of Fisher-Bingham distributions. The normalizing constant of the S1 has an alternative expression, as shown in the following.

Proposition 2**.**

For any $h>0$, let $\xi=(\frac{\nu(2\kappa_{0}+\kappa_{1})}{2(\kappa_{0}+h)},\frac{\kappa_{1}\sqrt{1-\nu^{2}}}{2h},0,\ldots,0)^{\top}\in\Re^{p}$ and let $\Psi$ be the $p\times p$ diagonal matrix with diagonal elements $(\kappa_{0}+h,h,\ldots,h)$. Moreover, let $g(r)$ $(r>0)$ be the probability density function of $R=Z^{\top}Z$, where $Z\sim N_{p}(\xi,\frac{1}{2}\Psi^{-1})$. Then the normalizing constant $a(\kappa_{0},\kappa_{1},\nu)$ of (5) is

[TABLE]

In Proposition 2, the function $g$ is the density of a linear combination of independent noncentral $\chi^{2}_{1}$ random variables. Following Kume and Wood (2005), we use saddle-point density approximations in the numerical computation of $g(1)$. First, note that the derivatives of the cumulant generating function, $K_{g}(t)=\log\int_{0}^{\infty}e^{tr}g(r)dr$, associated with the density $g$ have closed-form expressions. Denoting by $K_{g}^{(j)}(t)$ the $j$th derivative of $K_{g}(t)$, for $j=1,\ldots,4$, we get

[TABLE]

Let $\hat{t}$ be the unique solution in $(-\infty,h)$ of the saddle-point equation $K_{g}^{(1)}(t)=1$, which can be easily evaluated by using, e.g., a bisection method. Then a saddle-point density approximation of $g(1)$ is

[TABLE]

where $T={K_{g}^{(4)}(\hat{t})}/\{8(K_{g}^{(2)}(\hat{t}))^{2}\}-{5(K_{g}^{(3)}(\hat{t}))^{2}}/\{24{(K_{g}^{(2)}(\hat{t}))^{3}}\}$. In the following, we approximate the value of $a(\kappa_{0},\kappa_{1},\nu)$ by $\hat{a}(\kappa_{0},\kappa_{1},\nu)$ obtained by plugging (19) in place of $g(1)$ in (18).

We are now ready to describe our estimation procedure. Suppose $x_{1},\ldots,x_{n}$ is a sample from $\textrm{S1}(\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$. Given an initial value $\hat{\nu}^{(0)}\in(-1,1)$, we iterate between steps 1 and 2 below until convergence.

Step 1: Updating $\mu_{0}$ given all other parameters, of which we only need $\nu$. Suppose the inner product $\nu=\mu_{0}^{\top}\mu_{1}\in(-1,1)$ is fixed. Then maximizing the likelihood function with respect to $\mu_{0}$ is equivalent to minimizing $\frac{1}{n}\sum^{n}_{i=1}(\mu_{0}^{\top}x_{i}-\nu)^{2}$ subject to the constraint $\mu_{0}^{\top}\mu_{0}=1$. With a Lagrangian multiplier $\lambda$ using matrix notation, we solve

[TABLE]

where $\mathbb{X}$ is the $p\times n$ matrix whose $i$th column is $x_{i}$, yielding the necessary condition $S\mu_{0}-\nu\bar{x}-\lambda\mu_{0}=0,$ where $S=\mathbb{X}\mathbb{X}^{\top}/n$, $\bar{x}=\frac{1}{n}\sum^{n}_{i=1}x_{i}$. For a fixed Lagrangian multiplier $\lambda$, the solution is $\hat{\mu}_{0}=\nu(S-\lambda I_{p})^{-1}\bar{x}$, provided that $S$ is of full rank. The constraint $\mu_{0}^{\top}\mu_{0}=1$ makes us find a root $\lambda$ of $\nu^{2}\bar{x}^{\top}(S-\lambda I_{p})^{-2}\bar{x}-1$. The root $\hat{\lambda}$ is found by a bisection search in the range $[-\nu^{2}\bar{x}^{\top}\bar{x},\lambda_{S}]$, where $\lambda_{S}>0$ is the smallest eigenvalue of $S$ (Browne, 1967). The solution to (20) is then

[TABLE]

If $\nu=0$, then $\hat{\mu}_{0}$ is the eigenvector of $S$ corresponding to the smallest eigenvalue.

Step 2: Updating $(\mu_{1},\kappa_{0},\kappa_{1})$ given $\mu_{0}$. To facilitate the estimation of $\mu_{1}$, let $E(X)\in\Re^{p}$ be the usual expected value of the random vector $X$ in the ambient space $\Re^{p}$. We use the fact that $\mu_{1}$ is a linear combination of $\mu_{0}$ and $\gamma_{0}:=E(X)/\|E(X)\|$ (due to Lemma 3 in the Appendix) and reparameterize $\mu_{1}$ by an angle $\varphi$ and the direction ${\gamma^{*}}=P_{\mu_{0}}\gamma_{0}/\|P_{\mu_{0}}\gamma_{0}\|$ orthogonal to $\mu_{0}$, giving

[TABLE]

With $\gamma_{0}$ estimated by $\hat{\gamma}_{0}=\bar{x}/\|\bar{x}\|$, we now optimize for $\varphi\in[0,2\pi)$ together with $\kappa_{0},\kappa_{1}$, as follows. With $\mu_{0}$ and $\hat{\gamma}_{0}$ (thus $\hat{\gamma}^{*}$) given, the approximate negative log-likelihood with respect to $(\varphi,\kappa_{0},\kappa_{1})$ is

[TABLE]

Numerically minimizing (23) is much simpler than optimizing for $\mu_{1}$ with the nonlinear constraint $\|\mu_{1}\|=1$. We use a standard optimization package to obtain $\hat{\varphi},\hat{\kappa}_{0},\hat{\kappa}_{1}$ that minimizes (23). We get $\hat{\mu}_{1}$ by substituting $(\gamma^{*},\varphi)$ by $(\hat{\gamma}^{*},\hat{\varphi})$ in (22) and $\hat{\nu}=\cos(\hat{\varphi})$.

Let us now describe an extension of the above algorithm to the iMS1 model. Suppose $(x_{i(1)},\ldots,x_{i(K)})\in(S^{p-1})^{K}$ for $i=1,\ldots,n$ is a sample from an iMS1 model, where each marginal distribution is $\textrm{S1}(\mu_{0},\mu_{1(j)},\kappa_{0(j)},\kappa_{1(j)})$. While Step 2 above can be applied to update $\mu_{1(k)},\kappa_{0(k)},\kappa_{1(k)}$ given $\mu_{0}$, we modify Step 1 by replacing (20) with

[TABLE]

where the marginal $p\times n$ observation matrix $\mathbb{X}_{(j)}$ has the columns $x_{i(j)}$ $(i=1,\ldots,n)$. This is solved with the obvious analog to (21).

3.2 Estimation via profile likelihood for S2, iMS2 and MS2

The S2 model and its extensions have the convenient property that the horizontal components are independent of the vertical ones. To take advantage of this, suppose for now that $\mu_{0}$ is known. This allows to decompose an observation $x$ into two independent random variables $(s,y)$, which in turn leads to an easy estimation of the remaining parameters $\eta:=(\mu_{1},\kappa_{0},\kappa_{1})$. Thus our strategy of computing the m.l.e. proceeds in two nested steps. Let $\ell_{n}(\mu_{0},\eta)$ be the negative likelihood function given a sample $x_{1},\ldots,x_{n}$ from S2$(\mu_{1},\eta)$. In the outer step, we update $\mu_{0}$ to maximize a profile likelihood, i.e.,

[TABLE]

where evaluating

[TABLE]

for a fixed $\mu_{0}$ is the inner step. It is straightforward to see that the m.l.e. of $(\mu_{0},\eta)$ is given by $(\hat{\mu}_{0},\hat{\eta}_{\hat{\mu}_{0}})$.

In the following, we discuss in detail the inner step (25) of minimizing $\ell_{\mu_{0}}(\eta):=\ell_{n}(\mu_{0},\eta)$ for the iMS2 model (13) and for the MS2 model (16), while we resort to a standard optimization package for solving (24).

Independent multivariate S2 (iMS2). Suppose $(x_{i(1)},\ldots,x_{i(K)})\in(S^{p-1})^{K}$ for $i=1,\ldots,n$ is a sample from an iMS2, where each marginal distribution is S2$(\mu_{0},\mu_{1(j)},\kappa_{0(j)},\kappa_{1(j)})$. For a given $\mu_{0}$, the joint density can be written in terms of $(\mathbf{s}_{i},\boldsymbol{\phi}_{i})$ as done in (16), but with $\boldsymbol{\Lambda}=0$. Furthermore, with $\boldsymbol{\kappa}_{0}=(\kappa_{0(1)},\ldots,\kappa_{0(K)})^{\top}$, \mbox{\boldmath{\nu}}=(\nu_{(1)},\ldots,\nu_{(K)})^{\top}, we can write

[TABLE]

where $\circ$ denotes the element-wise product, and hence

[TABLE]

Note that the normalizing constant $b(\kappa_{0},\kappa_{1},\nu)$ satisfies

[TABLE]

where $\Phi(\cdot)$ is the standard normal distribution function. Finally, the negative log-likelihood function (given $\mu_{0}$) is

[TABLE]

where

[TABLE]

Therefore, the optimization for the inner step (25) is equivalent to simultaneously solving $2K$ subproblems.

Each of the $K$ subproblems of (3.2) is equivalent to obtaining the m.l.e. of a truncated normal distribution $\textrm{trN}(\nu_{(j)},(2\kappa_{0(j)})^{-1/2};(-1,1))$ based on the observations $s_{i(j)}$ ($i=1,\ldots,n$). Similarly, each of the $K$ subproblems of (29) amounts to obtaining the m.l.e. of a von Mises distribution with mean $\zeta_{(j)}$ and concentration $\kappa_{1(j)}$ from the sample $\phi_{i(j)}$ ($i=1,\ldots,n$). The m.l.e.s of the truncated normal are numerically computed, and we use the method of Banerjee et al. (2005) to obtain approximations of the m.l.e.s of the von Mises.

MS2. Under the general MS2 model (16) with a dependence structure on $\boldsymbol{\phi}_{i}$, a decomposition $\ell_{\mu_{0}}(\boldsymbol{\nu},\boldsymbol{\zeta},\boldsymbol{\kappa}_{0},\boldsymbol{\kappa}_{1},\boldsymbol{\Lambda})=\ell_{\mu_{0}}^{(1)}(\boldsymbol{\nu},\boldsymbol{\kappa}_{0})+\ell_{\mu_{0}}^{(2)}(\boldsymbol{\zeta},\boldsymbol{\kappa}_{1},\boldsymbol{\Lambda})$, similar to (27), is valid, where (29) is replaced by

[TABLE]

Minimizing (30) is equivalent to computing the m.l.e. of the multivariate von Mises distribution (Mardia et al., 2008). We either use maximum pseudo-likelihood estimate as discussed in Mardia et al. (2008) or moment estimates, yielding

[TABLE]

where $\bar{S}=(\bar{S}_{jk})$ and $\bar{S}_{jk}=\frac{1}{n}\sum_{i=1}^{n}\sin(\phi_{i(j)}-\hat{\zeta}_{(j)})\sin(\phi_{i(k)}-\hat{\zeta}_{(k)})$ for $j,k=1,\ldots,K$. These estimates coincide with the m.l.e.s when $K=2$. For larger $K>3$, the accuracy of the moment estimates deteriorates, but evaluating m.l.e.s or a maximum pseudo-likelihood estimator becomes computationally highly expensive.

4 Testing hypotheses

It is of interest to infer on the parameters of our models. In this section, we describe a large-sample testing procedure for several hypotheses of interest.

Our testing procedure is based on the likelihood ratio statistic, with effort devoted to an identification of the restricted parameter space $\Theta_{0}$, for each hypothesis and computing the maximized likelihood under $\Theta_{0}$. Recall that the parameter space for the iMS1 and iMS2 models is given by $\Theta_{\rm ind}=\mathbb{S}^{p-1}\times(\mathbb{S}^{p-1})^{K}\times(\Re_{+})^{K}\times(\Re_{+})^{K}$ for \theta_{\rm ind}=(\mu_{0},\mbox{\boldmath{\mu}}_{1},\mbox{\boldmath{\kappa}}_{0},\mbox{\boldmath{\kappa}}_{1}). For the more general GMS2 model including associations, we have $\Theta_{\rm GMS2}=\Theta_{\rm ind}\times(\Re^{(p-1)^{2}})^{K(K-1)/2}$ for $\theta_{\rm GMS2}=(\theta_{\rm ind},\mathbf{B})$. In the following, we describe our testing procedure using the MS2 distribution in dimension $p=3$, whose parameter space is $\Theta=(\theta_{\rm ind}\times(\Re)^{K(K-1)/2})$ for \theta=(\theta_{\rm ind},\mbox{\boldmath{\Lambda}}). For some $\Theta_{0}$ that dictates a null hypothesis $H_{0}$ and satisfies $\Theta_{0}\subset\Theta$, we denote the maximized log-likelihood under $\Theta_{0}$ by $\mathcal{L}_{0}$, and the maximized log-likelihood under $\Theta$ by $\mathcal{L}_{1}$. It is well-known that for this nested model, for large sample size $n$, $W_{n}:=-2(\mathcal{L}_{0}-\mathcal{L}_{1})$ follows approximately a chi-square distribution with $q_{1}-q_{2}$ degrees of freedom, where $q_{1}=(p-1)(K+1)+2K+K(K-1)/2$ and $q_{2}$ are the dimensions of $\Theta$ and $\Theta_{0}$, respectively. Once $W_{n}$ is computed, as usual, our test rejects the null hypothesis for large enough values of $W_{n}$.

We are interested in the following null hypotheses, with the alternative being the full MS2 distribution. Testing the first hypothesis gives a test of association among directional vectors, while the latter two provide tests regarding the parameters of the underlying small-sphere \mathcal{C}(\mbox{\boldmath{\mu}}_{0},\mbox{\boldmath{\nu}}). In all three cases, the alternative is $H_{1}:\theta\in\Theta\setminus\Theta_{0}$.

1. Test of association. $H_{0}$: $\boldsymbol{\Lambda}=\textbf{0}$, i.e., $\theta\in\Theta_{0}=\mathbb{S}^{p-1}\times(\mathbb{S}^{p-1})^{K}\times(\Re_{+})^{K}\times(\Re_{+})^{K}\times\{\textbf{0}\}$. Under $H_{0}$, the model degenerates to the iMS2 and there is no horizontal dependence.

2. Test of axis. $H_{0}$: $\mu_{0}=\mu_{0}^{*}$, i.e., $\theta\in\Theta_{0}=\{\mu_{0}^{*}\}\times(\mathbb{S}^{p-1})^{K}\times(\Re_{+})^{K}\times(\Re_{+})^{K}\times(\Re)^{K(K-1)/2}$. This is to test whether a predetermined axis $\mu_{0}^{*}$ of the small sphere is acceptable.

3. Test of great-sphere. $H_{0}$: \mbox{\boldmath{\nu}}=0, i.e., $\theta\in\Theta_{0}\simeq\mathbb{S}^{p-1}\times(\mathbb{S}^{p-2})^{K}\times(\Re_{+})^{K}\times(\Re_{+})^{K}\times(\Re)^{K(K-1)/2}$. ($A\simeq B$ means that $A$ and $B$ are diffeomorphic.) This is to test whether the underlying spheres are great spheres with radius 1.

While the test of association (Hypothesis 1) is only available under the MS2 model ($p=3$), Hypotheses 2 and 3 can also be tested using S1, S2, iMS1 or iMS2 models in any dimension $p\geq 3$. Moreover, to validate the use of small-sphere distributions, in any dimension $p\geq 3$, the following hypotheses can be tested. For simplicity, assume for now the S1 model with $\theta=(\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})\in\Theta=(\mathbb{S}^{p-1})^{2}\times(\Re_{+})^{2}$.

4. Test for von Mises-Fisher distribution. $H_{0}$: $\kappa_{0}=0$, i.e., $\theta\in\Theta_{0}\simeq\mathbb{S}^{p-1}\times\Re_{+}$. Under $H_{0}$, there is no “small-circle feature.”

5. Test for Bingham-Mardia distribution. $H_{0}$: $\kappa_{1}=0$, i.e., $\theta\in\Theta_{0}\simeq\mathbb{S}^{p-1}\times\Re_{+}$. Under $H_{0}$, there is no unique mode.

Hypotheses 4 and 5 can also be tested under the S2, iMS1, iMS2 and MS2 models. Additional care is needed when using S2, iMS2, or MS2 for Hypothesis 4, as the null distribution is the von Mises-Fisher on \mathcal{C}(\mbox{\boldmath{\mu}}_{0},\mbox{\boldmath{\nu}}) (not on $\mathbb{S}^{p-1}$).

For each hypothesis, computing the test statistic $W_{n}$ requires to maximize the likelihood on $\Theta_{0}$ (or to compute $\mathcal{L}_{0}$). This is easily achieved by modifying the iterative algorithms in Section 3. For example, for the test of association, computing $\mathcal{L}_{0}$ and $\mathcal{L}_{1}$ amounts to obtaining the m.l.e.s under the iMS2 and MS2 models, respectively; for Hypothesis 2 (test of axis), where $\mu_{0}$ is given, one only needs to solve (25) once. Other cases of restricted m.l.e.s can be easily obtained. In the online supplementary material, we confirm that the test statistic $W_{n}$ using our algorithms under the null hypotheses above are empirically nearly chi-square distributed for sample size $n=30$. In Section 5.3 and in the online supplementary material, empirical powers of the proposed test procedures are reported for several important alternatives.

5 Numerical studies

We demonstrate the performances of small-circle fitting in Section 5.1, the ability of the MS2 of modeling the horizontal dependence in Section 5.2 and a testing procedure to prevent overfitting in Section 5.3.

5.1 Estimation of small-circles

The performance of our estimators in fitting the underlying small-spheres $\mathcal{C}(\mu_{0},\nu)$ is numerically compared with those of competing estimators obtained from assuming the Bingham-Mardia (BM) distribution and the least-square estimates of Schulz et al. (2015). The BM distribution has originally been defined only for data on $\mathbb{S}^{2}$, but we use a natural extension given by a special case of the iMS1. Thus, “BM estimates” refer to the estimates of the iMS1 model with the restriction $\boldsymbol{\kappa}_{1}=0$. The estimates of Schulz et al. are obtained by minimizing the sum of squared angular distances from observations to $\mathcal{C}(\hat{\mu}_{0},\hat{\nu})$, which will be referred to as a “least-squares (LS)” method.

We first consider four univariate S2 models to simulate data concentrated on a small circle. The directional parameters $(\mu_{0},\mu_{1})$ are set to satisfy $\nu=0.5$. We use $(\kappa_{0},\kappa_{1})=(10,1),(100,1),(100,10)$ to represent various data situations. Random samples from these three settings are shown in Fig. 4(a)–(c). We also consider the BM model as a special case of the S2 distributions (by setting $\kappa_{1}=0$); a sample from the BM distribution is shown in Fig. 4(d).

The small-circle estimation performances of the S1, S2, BM, and LS estimates are measured by an angular product error (in degrees), defined as

[TABLE]

Table 2 displays the means and standard deviations of ${L}\left((\mu_{0},\nu),(\hat{\mu}_{0},\hat{\nu})\right)$ from 100 repetitions for each of the four methods, fitted to random samples of size 50 from each of the settings, labeled (a)–(d). The S2 estimates performed best for all non-trivial settings: (a), (b) and (c). Estimators under Settings (a) and (c) show higher standard errors. This is expected and due to either the large vertical dispersion or the smaller horizontal dispersion as visualized in Fig. 4(a) and (c). Even when the sophisticated S1 and S2 models are not needed in Setting (d), the S1 and S2 estimators perform virtually as good as the BM estimator does.

Next, to show the performance of our multivariate models, we consider six bivariate MS2 models. The directional parameters (\mu_{0},\mbox{\boldmath{\mu}}_{1}) were set to satisfy \mbox{\boldmath{\nu}}=(0.5,-.3), and the concentration parameters were chosen to mimic the concentrations of the univariate models, described above. For Cases (a)–(c), we set $(\kappa_{0j},\kappa_{1j},\lambda_{12})=(10,1,0),(100,1,0),(100,10,0)$, for $j=1,2$, so that the models are indeed the iMS2. For the later three cases (d)–(f), we set $(\kappa_{0j},\kappa_{1j},\lambda_{12})=(10,2,1.5),(100,2,1.5),(100,20,15)$, $j=1,2$, to make their vertical and horizontal dispersions be similar to the iMS2 counterparts. By setting $\lambda_{12}>0$, the random bivariate directions are positively associated. (Examples of random samples from these settings can be found in Fig. A6 in the online supplementary material.) The small-circle estimation performance of the iMS1, iMS2, MS2, BM and LS estimates is measured by the canonical multivariate extension of the angular product error (32). Table 3 collects the means and standard deviations of the angular product errors from 100 repetitions with the sample size $n=50$. Overall, the three proposed models (iMS1, iMS2, and MS2) show better or competitive performances in the axis and radii estimation. In particular, when directions are clearly concentrated on small-circles and are horizontally dependent, i.e., in Settings (e) and (f), the MS2 estimates shows better performances than others.

We check robustness against model misspecification of the estimators by simulating data from a more general signal-plus-noise model (neither S1 nor S2). The performances of small-circle fitting of the proposed methods are comparable to that of the least-square estimator. Relevant simulation results and a detailed discussion can be found in the online supplementary material.

5.2 Estimation of horizontal dependence

The ability of the MS2 to model the horizontal dependence is an important feature of the proposed distributions. Here, we empirically confirm that the MS2 estimates provide accurate measures of horizontal dependence, using Cases (c) and (f) in Section 5.1. For sample sizes $n=50$ and $200$, the concentration and association parameters were estimated under the assumption of MS2 (or iMS2), and Table 4 summarizes the estimation accuracy. In all cases, the MS2 model provides precise estimations of the horizontal dispersion and dependence; as the sample size increases, the mean squared error decreases. For Case (c), the underlying model is exactly iMS2, so the iMS2 estimates have smaller mean squared errors than the MS2 estimates. However, for Case (f), we notice that the iMS2 estimates of $\boldsymbol{\kappa}_{1}=(\kappa_{11},\kappa_{12})$ become inferior. In fact, in case of existing horizontal dependence, i.e., when $\lambda_{12}\neq 0$, the concentration parameters \mbox{\boldmath{\kappa}}_{1} in the misspecified iMS2 model do not correctly represent the concentrations as correctly represented by the MS2 model. This is so because the marginal distribution of $\phi_{j}$, $j=1,2$, in (16) is not a von Mises distribution (Shing et al., 2002; Mardia et al., 2008).

5.3 Detecting overfitting in an isotropic case

When the data do not exhibit a strong tendency of a small-circle feature, the S1 and S2 distributions may overfit the data. For example, to a random sample from an isotropic vMF distribution as shown in Fig. 5(a), the S1 or S2 model fits an unnecessary small-circle $\mathcal{C}(\hat{\mu}_{0},\hat{\nu})$. Indeed, a small-circle fit was observed in 83% of simulations of fitting the S1 model. Using the BM model or the LS results in a similar overfitting, where very small circles are erroneously fitted for 100% and 68% of the simulations, for the BM and LS, respectively.

This problem of overfitting has been known for a while and discussed in the context of dimension reduction of directional data. In particular, Jung et al. (2011, 2012) and Eltzner et al. (2015) investigated the overfitting phenomenon for the least-square estimates and proposed some ad-hoc methods for adjustment. To prevent the overfitting, we point out that the testing procedure in Section 4 for the detection of isotropic distributions (Hypothesis 4) works well. To confirm this, we evaluated the empirical power of the test at the significance level $\alpha=0.05$ for several alternatives. The power increases sharply as the distributions become more anisotropic; under the alternative distributions depicted in Fig. 5(b)–(d), the empirical powers are respectively $\hat{\beta}=0.435,1$ and $1$, evaluated from 200 repetitions.

6 Analysis of s-rep data

In this section, an application of the proposed distributions and test procedures to s-rep data is discussed.

6.1 Modeling rotationally-deformed ellipsoids via s-reps

Skeletal representations (s-reps) have been useful in mathematical modeling of human anatomical objects (Siddiqi and Pizer, 2008). Roughly, an s-rep model for a 3-dimensional object consists of locations of a skeletal mesh (inside of the object) and spoke vectors (directions and lengths), connecting the skeletal mesh with the boundary of the object; examples are shown in the top left panel of Fig. 6. When the object is “rotationally deformed”, Schulz et al. (2015) have shown that the directional vectors of an s-rep model approximately trace a set of concentric small-circles on $\mathbb{S}^{2}$, as shown in the top panels of the figure. Such rotational deformations (e.g., rotation, bending and twisting) of human anatomical objects have been observed in between and within shape variations of hippocampi and prostates (Joshi et al., 2002; Jung et al., 2011; Pizer et al., 2013). We demonstrate a use of the MS2 distribution in modeling (and fitting) a population of such objects via s-reps. Note that the sample space of an s-rep with $K$ spokes is $(\mathbb{S}^{2})^{K}\times\Re_{+}^{K}\times(\Re^{3})^{K}$ (for direction, length, and location). In this work, we choose to analyze the spoke directions in $(\mathbb{S}^{2})^{K}$ only, leaving a full-on analysis, accommodating the lengths and locations, to future work.

6.2 Data preparation

For our purpose of validating the use of the proposed distributions, we use an s-rep data set, fitted from $30$ deformed ellipsoids; two samples from this data set are shown in Fig. 6(a). This data set was previously used in Schulz et al. (2015) as a simple experimental representation of real human organs. The data set was generated by “physically bending” a template ellipsoid about an axis $\mu^{*}_{0}=(0,1,0)$ by random angles drawn from a normal distribution with standard deviation 0.4 (radians). Each deformed ellipsoid is then recorded as a 3-dimensional binary image. To mimic the procedure of fitting s-reps from, for example, medical resonance imaging of a real patient, s-reps with $74$ spokes were fitted to these binary images. (See Pizer et al. (2013) for details of the s-rep fitting.) As a preprocessing, we chose $K=58$ spoke vectors, excluding the vectors with very small total variation.

6.3 Inference on the bending axis

Fitting the iMS1 distribution, we obtained axis estimate $\hat{\mu}_{0}^{(\mathrm{iMS1})}=(0.007,1.000,-0.008)$ (rounded to three decimal digits). Similarly, from the MS2 fitting, $\hat{\mu}_{0}^{(\mathrm{MS2})}=(0.006,1.000,0.006)$. These estimates are virtually the same, only $0.6$ degrees away from the ground truth $\mu_{0}^{*}$. Estimates of the concentric small-circles $\mathcal{C}(\hat{\mu}_{0}^{(\mathrm{MS2})},\hat{\nu}_{j})$ for four choices of $j$ (the spoke index) are also shown in the top right panel of Fig. 6 in which $\hat{\mu}_{0}^{(\mathrm{MS2})}$ and $\mu_{0}^{*}$ are also shown. The test of axis (Hypothesis 2) discussed in Section 4 is applied to test $H_{0}:\,\mu_{0}=\mu_{0}^{*}$ under the iMS2 model. With the test statistic assuming the value $W_{n}=2.69$, compared to 5.99 (the 95th percentile of $\chi_{2}^{2}$), and the corresponding p-value of 0.26, we cannot reject that the s-rep spokes are rotated about the bending axis $\mu_{0}^{*}$.

6.4 Inference on horizontal dependence

An advantage of modeling the s-rep spoke directions by the MS2 distribution is the ability of perceiving and modeling the horizontal dependence among directions. As an exploratory step, we have collected the estimated correlation coefficients, computed from the approximate precision matrix $\widehat{\Sigma}^{-1}$, whose elements are \hat{}\mbox{\boldmath{\kappa}}_{1} and \widehat{}\mbox{\boldmath{\Lambda}}; see (17). A histogram of $K(K-1)/2$ estimated correlation coefficients is plotted in the bottom left panel of Fig. 6. Notably, pairs of spoke vectors from the same side (e.g. two spoke vectors in the “left side” of the ellipsoids in Fig. 6) exhibit strong positive correlations, while those from the opposite sides exhibit strong negative correlations. The horizontal dependence is in fact apparent by the way data were generated (simultaneously bending all the spoke directions).

For large enough sample sizes, we could use the test of association discussed in Section 4 for testing $H_{0}:\,\boldsymbol{\Lambda}=\textbf{0}$. Unfortunately, due to our small sample size, $n=30$, and the large number of parameters tested, 1653 ($=K(K-1)/2$), this is infeasible. Coping with this high-dimension, low-sample-size situation is beyond the scope of current paper, and we resort to choose only two spoke directions to test the dependence, but to repeat the testing for many different pairs of total $K=58$ spokes. For each pair of spokes, the likelihood-ratio test produces a p-value for the pair. Investigating the empirical distribution of these p-values can provide a rough estimate of the power. In Fig. 6(d), it can be seen that, at the significance level $0.05$, the MS2 test of dependence is indeed powerful, with a rejection rate of 97%.

To provide some context to this rate, the MS2 test was compared with other natural choices of tests. We applied two methods that were previously used for s-rep data analysis: the composite principal nested spheres (CPNS), discussed in Pizer et al. (2013), and the least-square (concentric) small-circle fitting method of Schulz et al. (2015).

The CPNS-test is built as follows. First, the least-square small-circle is fitted to each marginal directions on $\mathbb{S}^{2}$. With an understanding that the axis of the fitted small-circle points to the north pole, the observations (say, $x_{i(k)}$ from the $i$th sample, $k$th spoke) are represented in spherical coordinates ($\theta_{i(k)},\phi_{i(k)}$). For the purpose of testing “horizontal associations”, we only keep the longitudinal coordinates $\theta_{i(k)}$. For any given pair ($k,\kappa$), Fisher’s z-transformation is used to obtain the p-value in testing whether the correlation coefficient between $\theta_{i(k)}$ and $\theta_{i(\kappa)}$ is zero. We refer to this test procedure by a CPNS test.

An LS test procedure is defined similarly to the CPNS test, except that the first step of fitting individual small-circles is replaced by fitting concentric small-circles.

These two tests were also conducted for the same combinations of spoke directions, and the empirical distributions of respective p-values are also plotted in Fig. 6. These alternative tests appear to be too conservative, with rejection rates 11% for the LS test, and 13.5% for the CPNS test (at level 0.05). Heuristically, the higher power of the MS2 tests is due to the superior fitting of the MS2 distribution. In particular, the “horizontal angles” predicted from the MS2 tend to be linearly associated, while those from the least-squares fit tend to be arbitrary. We refer to the online supplementary material for more simulation results. All in all, using the MS2 distribution shows a clear advantage in modeling and testing the horizontal dependence of multivariate directions.

7 Human knee gait analysis

In biomechanical gait analysis, accurately modeling human knee motion during normal walking has a potential to differentiate diseased subjects from normal subjects. In particular, the axis of bending (of the lower leg toward the upper leg) is believed to be a key feature in the discrimination among the diseased and normal subjects (Pierrynowski et al., 2010). As a step towards the development of statistical tests for a “two-group axes difference,” in this section we employ the proposed distributional families in modeling the bending motion of the knee.

The raw data set we use is obtained from a healthy volunteer and it is a time series of coordinates of markers planted at the volunteer’s leg, recorded for 16 gait cycles. For each time point, the directional vectors on $(\mathbb{S}^{2})^{5}$ were computed to be the unit vectors between reference markers, as done in Schulz et al. (2015). These directional vectors are horizontally dependent of each other (as evidenced in Fig. 7), which suggests that we may fit the MS2 distribution.

The first panel of Fig. 7 illustrates the result of MS2 fit to the all data points. There, we superimpose the fitted concentric circles to the observed directional vectors, including their estimated axis, together with a hypothesized dominant bending axis $\mu_{0}^{*}=(0,1,0)^{\top}$, the left-right axis of the subject. The MS2 model seems to fit well with high estimated horizontal correlation coefficients. We, however, identify a strong evidence against using a single MS2. Specifically, as shown in Fig. 7 some directional vectors exhibit higher variations for a subset of time points.

In fact, the data consist of many inhomogeneous periods of the gait cycle. We focus on the “swing” and “stance” periods, and separately analyze subsampled data from each period. The MS2 model fits well for the swing period data (see Fig. 7(b)), and the estimated axis $\hat{\mu}_{0}^{(\mbox{sw})}=(0.013,1.000,0.005)^{\top}$ is only 0.8 degrees away from the hypothesized axis, $\mu_{0}^{*}$. As expected, our likelihood ratio test procedure does not reject the null hypothesis $H_{0}:\mu_{0}=\mu_{0}^{*}$, with p-value 0.16. For the stance period data, excluding the highly-irregular directions shown as dark blue points in Fig. 7(c), the MS2 model also fits well, and we confirm that the axis of bending for this period differs from $\mu_{0}^{*}$, with p-value less than $10^{-5}$. The estimated axis for the stance period is $\hat{\mu}_{0}^{(\mbox{st})}=(0.11,0.994,0.006)^{\top}$.

While the MS2 distribution was useful in making inference on the bending axis of partial knee motions, future work for this type of data lies in the development of a two-sample axis difference test.

Acknowledgements

Stephan Huckemann gratefully acknowledges funding by the Niedersachsen Vorab of the Volkswagen Foundation and DFG HU 1575/4.

Supporting information

Additional information for this article is available online. These include graphics of S2 densities (Appendix A1), powers of the test (Appendix A2), and additional simulation results (Appendix A3 and A4).

Address

Department of Statistics, University of Pittsburgh, 1806 Wesley W. Posvar Hall, 230 Bouquet Street, Pittsburgh, PA 15260. sungkyu@pitt.edu

Appendix

We provide a technical lemma, referenced in Section 3.1, and proofs of Propositions 1 and 2.

Lemma 3**.**

If $X\sim$S1($\mu_{0},\mu_{1},\kappa_{0},\kappa_{1}$), then $E(X)$ is a linear combination of $\mu_{0}$ and $\mu_{1}$.

Proof of Lemma 3.

Suppose that for some $a,b,c\in\Re$, $v\in S^{p-1}$, $E(X)=a\mu_{0}+b\mu_{1}+cv$. Then choose a $B\in O(p)$ such that $B\mu_{0}=\mu_{0}$, $B\mu_{1}=\mu_{1}$ but $Bv\neq v$. By Proposition 1(i), $BX\sim S1(\mu_{0},\mu_{1},\kappa_{0},\kappa_{1})$. Thus $E(X)=E(BX)$, which in turn leads to $a\mu_{0}+b\mu_{1}+cv=a\mu_{0}+b\mu_{1}+cBv$, which is true only if $c=0$. This gives the result. ∎

Proof of Proposition 1.

Assertions (i) and (iii) are trivial. (ii): Assume that $X\sim BX$. With the representation (10) this is equivalent with

[TABLE]

In consequence

[TABLE]

Thus the case of $B\mu_{0}=\mu_{0}$ yields at once $B\mu_{1}=\mu_{1}$. In contrast, the case $B\mu_{0}=-\mu_{0}$ yields that $\|B\mu_{1}\|^{2}=\|4\nu\frac{\kappa_{0}}{\kappa_{1}}\mu_{0}+\mu_{1}\|^{2}=1+8\nu^{2}\kappa_{0}/\kappa_{1}(1+2\kappa_{0}/\kappa_{1})>1$, which cannot be if $B$ is orthogonal.

Further, if $Y\sim BY$, since the exponent is again the sum of an even function in $x$ and an odd function, we have again that $B\mu_{0}=\pm\mu_{0}$, yielding $\|P_{\mu_{0}}x\|=\|P_{\mu_{0}}Bx\|$, due to orthogonality, $B$ preserves the space orthognal to $\mu_{0}$, as well as

[TABLE]

As before, $B\mu_{0}=\mu_{0}$ yields that $B\mu_{1}=\mu_{1}$ and $B\mu_{0}=-\mu_{0}$ would give that

[TABLE]

In particular, this would yield that $B\mu_{1}-\mu_{1}$ is a multiple of $\mu_{0}$, the factor, however, is not constant in $x$, a contradiction. ∎

Proof of Proposition 2.

For a given $h>0$, let $\gamma=2\kappa_{0}\nu\mu_{0}+\kappa_{1}\mu_{1}$ and $A_{h}=\kappa_{0}\mu_{0}\mu_{0}^{\top}+hI_{p}$. Then it is easy to see that the S1 density (5) can be expressed as the Fisher-Bingham form (10):

[TABLE]

where $\alpha(\gamma,A_{h})$ satisfies

[TABLE]

For the purpose of evaluating the value of $a(\kappa_{0},\kappa_{1},\nu)$, or equivalently $\alpha(\gamma,A_{h})$ for the given value of $h$, one can assume without losing generality that $\mu_{0}=(1,0,\ldots,0)^{\top}$ and $\mu_{1}=(\nu,\sqrt{1-\nu^{2}},0,\ldots,0)$, so that $\gamma=(\nu(2\kappa_{0}+\kappa_{1}),\kappa_{1}\sqrt{1-\nu^{2}},0,\ldots,0)^{\top}$, and the vector of diagonal values of $A_{h}$ are $\lambda:=(2(\kappa_{0}+h),h,\ldots,h)$. The $j$th element of $\xi$, in the statement of proposition, is then given by $\mu_{j}:=\gamma_{j}/2\lambda_{j}$. With these notations, Proposition 1 of Kume and Wood (2005) gives

[TABLE]

Hence, by (33), we have (18). ∎

See pages 1-16 of SSD-supp-2017-05-27.pdf