10.2.1 The density

We consider the case where we have a probability distribution on the pre-shape sphere S^{m(k − 1) − 1}, corresponding to k points in m dimensions. For the m = 2 dimensional case and using complex notation we have seen in Section 4.3 that , where is the unit complex sphere in k − 1 complex dimensions. In our particular application to shape analysis consider the k original landmarks in m = 2 dimensions with complex coordinates written as z^o (k × 1 complex vector). Pre-multiply by the (k − 1) × k Helmert submatrix H of Equation (2.10) to give the k − 1 Helmertized complex landmarks z_H [a (k − 1)-vector]. Rescaling by ||z_H|| it follows that the pre-shape is:

Definition 10.1 The complex Bingham distribution on , denoted , has probability density function (p.d.f.)

where z* represents the complex conjugate of the transpose of z, the matrix A is (k − 1) × (k − 1) Hermitian (i.e. A = A*), and c(A) is the normalizing constant.

The complex Bingham distribution is analogous to the real Bingham distribution (e.g. see Mardia et al. 1979, p. 433) and in fact, it is a particular case (see Section 10.2.5). The distribution has the property that

and is thus invariant under rotations of the pre-shape z. So, if an object is rotated, then it has the same density and will contribute in an identical way to inference as the same object in the original rotation. This property therefore makes the distribution appropriate for shape analysis (location and scale were previously removed because z is on the pre-shape sphere). The complex Bingham distribution provides a very elegant framework for the analysis of 2D shape data. The distribution was proposed by Kent (1994) and the properties that we describe are from his paper.

Since z*z = 1 for , we can see that the parameter matrices A and A + αI define the same complex Bingham distribution, with c(A + αI) = c(A)exp α, where α is a complex number. It is convenient to remove this non-identifiability by setting λ_max (A) = 0, where λ_max (A) denotes the largest eigenvalue of A.

10.2.2 Relation to the complex normal distribution

In order to understand the covariance structure of a complex distribution we describe the complex multivariate normal distribution. If z_j = x_j + iy_j have a joint complex normal distribution with means ξ_j = μ_j + iν_j, j = 1, …, p, and a p × p Hermitian covariance matrix Σ = Σ₁ + iΣ₂, then writing x = (x₁, …, x_p, y₁, …, y_p)^T and μ = (μ₁, …, μ_p, ν₁, …, ν_p)^T we have

(10.3)

where Σ₂ = −Σ^T₂ is skew-symmetric and Σ₁ is symmetric positive definite. In particular, var(x_j) = var(y_j) and cov(x_j, y_j) = 0, and so for each j the covariance structure is isotropic. Writing z = (z₁, …, z_p)^T and ξ = (ξ₁, …, ξ_p)^T the density of the complex normal distribution is:

(10.4)

and we write ; see, for example, Goodman (1963) and Andersen et al. (1995).

The complex Bingham distribution can be obtained through conditioning a zero mean complex multivariate normal distribution to have unit norm. In particular, if then,

So, an interpretation of the Hermitian matrix − A is that of the inverse covariance matrix of a zero mean complex normal random variable, which is conditioned to have unit norm to give the complex Bingham distribution. Hence, the complex Bingham distribution is an example of a conditional approach.

10.2.3 Relation to real Bingham distribution

The (k − 2) dimensional complex Bingham distribution can be regarded as a special case of a (2k − 3) dimensional real Bingham distribution as in the following. If the jth element of z is (z)_j = x_j + iy_j, define a (2k − 2) dimensional vector u = (x₁, y₁, …, x_{k − 1}, y_{k − 1})^T, say, by splitting each complex number into its real and imaginary parts. Also, if A = (a_hj) has entries a_hj = α_hjexp (iθ_hj) with θ_jh = −θ_hj, − π < θ_jh ≤ π, define a (2k − 2) × (2k − 2) matrix B made up of (k − 1)² blocks of size (2 × 2) given by:

Then z*Az = u^TBu so that a complex Bingham distribution for z is equivalent to a real Bingham distribution for u (Kent 1994).

10.2.4 The normalizing constant

Result 10.3 (Kent 1994) The normalizing constant of the complex Bingham distribution is:

(10.5)

where λ₁ < λ₂ < ⋅⋅⋅ < λ_{k − 1} = 0 denote the eigenvalues of A. Note that c(A) = c(Λ) depends only on the eigenvalues of A, with Λ = diag(λ₁, …, λ_{k − 1}).

The result and proof were given by Kent (1994). If equalities exist between some λ_js then c(A) can be integrated using repeated applications of L’Hôpital’s rule. For example, for λ_{k − 2} = λ_{k − 1} but all other λ_js distinct, we get (Bingham et al. 1992):

A particularly useful and simple case is the complex Watson distribution when there are just two distinct eigenvalues, as described in Section 10.3.

10.2.5 Properties

We give various properties of the complex Bingham distribution. Let γ₁, …, γ_{k − 1} denote the standardized eigenvectors of A so that γ*_jγ_j = 1, γ*_iγ_j = 0, i ≠ j, and Aγ_j = λ_jγ_j, j = 1, …, k − 1. Each γ_j is defined only up to rotation by a unit complex scalar. Provided that λ_{k − 2} < 0, z = γ_{k − 1} maximizes the density and γ_{k − 1} is unique up to a scalar rotation exp (iθ)γ_{k − 1}. Furthermore, if λ₁, …, λ_{k − 2} are far below 0, the distribution becomes highly concentrated about this modal axis.

10.2.5.1 High concentrations

To model high concentration, replace A by κA transforming to Kent’s polar coordinates and integrating out θ₁, …, θ_{k − 1}

As κ → ∞, u_j = κs_j, j = 1, …, k − 2, tend to independent exponential random variables with densities

Thus for large κ,

(10.6)

We assume that there is a unique largest eigenvalue λ_{k − 1} = 0 > λ_{k − 2}. Another way to look at the high concentration about the axis γ_{k − 1} is to rotate z to γ_{k − 1} and then to project z onto the tangent plane at γ_{k − 1}. Given ,

and it can be proved that θ and v are independently distributed with θ uniform on [0, 2π). Writing A = ∑^{k − 1}_{j = 1}λ_jγ_jγ*_j = ∑^{k − 2}_{j = 1}λ_jγ_jγ*_j, we have

and therefore v is asymptotically complex multivariate normal in k − 1 complex dimensions with zero mean and covariance matrix

where the last expression is the Moore–Penrose generalized inverse of − A (see Definition 9.1). Here Σ is singular because v lies in the tangent plane γ*_{k − 1}v = 0. Thus z is approximately distributed as so that Σ is determined by the eigenvectors orthogonal to γ_{k − 1}.

The asymptotic distribution shows that the real and imaginary parts of v_i are jointly distributed as N₂(0, I₂σ²_i). Hence the complex Bingham distribution imposes an isotropic distribution on the marginal distributions. However, it does allow some intercorrelation between landmarks, that is that of a complex normal covariance matrix as seen in Equation (10.3).

10.2.6 Inference

Let z₁, …, z_n be a random sample from a population modelled by a complex Bingham distribution, n ≥ k − 1. Set

to be the (k − 1) × (k − 1) complex sum of squares and products matrix. Suppose that the eigenvalues of S are positive and distinct, 0 < l₁ < ⋅⋅⋅ < l_{k − 1} and let g₁, …, g_{k − 1} denote the corresponding eigenvectors. Note that ∑_jl_j = n. Maximum likelihood estimation is carried out as for the real Bingham distribution.

Result 10.4 (Kent, 1994) Under the complex Bingham distribution the MLE of γ₁, …, γ_{k − 1} and Λ = diag(λ₁, …, λ_{k − 1}) are given by:

and the solution to

Under high concentrations

Proof: The log-likelihood for the data reduces to:

(10.7)

Holding λ₁ < … < λ_{k − 2} < λ_{k − 1} = 0 constant, it can be seen that

and so

The MLEs of the eigenvalues are found by solving

Under high concentration, from Equation (10.6),

giving ,    j = 1, …, k − 2.

The dominant eigenvector can be regarded as the average axis of the data: an estimate of modal shape. Hence, the average axis from the complex Bingham MLE is the same as the full Procrustes estimator of mean shape for m = 2 dimensional data, which was also given by the dominant eigenvector of S as shown in Result 8.2.

Example 10.1 Consider the mouse vertebral data described in Section 1.4.1. There are k = 6 landmarks in m = 2 dimensions here. Using the complex Bingham distribution the MLE of the population modal shape is:

corresponding to the eigenvalue l_{k − 1} = 22.9. Since l_{k − 1} ≈ n the data are highly concentrated here. In Figure 8.5 we saw a plot of the MLE of the modal shape (which is the full Procrustes mean) as a centred pre-shape.

In order to examine the structure of variability we could examine g₁, …, g_{k − 2} although the restrictions of a complex covariance structure are in force. In particular, we necessarily have isotropic perturbations at each landmark. So, a practical method for inference would be to examine the shape variability in the tangent space; details of this example were given in Example 8.2.

10.2.7 Approximations and computation

In order to carry out more accurate estimation of the complex Bingham parameters Kume and Wood (2007) study the derivatives of the normalizing constant of the Bingham distribution, and Kume and Wood (2005) consider highly accurate saddlepoint approximations for the Bingham normalizing constant (as well as the Fisher–Bingham). Their methodology enables more accurate maximum likelihood estimation for a whole range of parameter values, including when the data are quite dispersed.

10.2.8 Relationship with the Fisher–von Mises distribution

Result 10.5 For triangles (k = 3) modelling with the complex Bingham in pre-shape space is equivalent to using the Fisher–von Mises distribution on the spherical shape space .

Proof: If k = 3, then there is just one non-zero eigenvalue of A. Hence we can re-parameterize A as:

So, we have

where ρ( · ) is the Riemannian distance. We can obtain the density in shape space by using Kent’s polar shape coordinates and integrating over the independent uniform θ₂ component (Mardia 1996b). So, the density in shape space is proportional to:

where H^Tz is the centred icon corresponding to z. Since the shape space can be regarded as a sphere with great circle metric ρ, it follows that, using the direction vector l = (l_x, l_y, l_z)^T, the probability element of l is:

(10.8)

where ν = (ν_x, ν_y, ν_z)^T is the modal direction and dS is the probability element. The density (10.8) is the density of the Fisher–von Mises distribution on .

Mardia (1989b) provided motivation to use the Fisher–von Mises distribution for shape analysis of triangles.

Example 10.2 Consider the microfossil data described in Section 1.4.16. Using the complex Bingham distribution on the pre-shape sphere will be equivalent to using the Fisher–von Mises distribution on . Let l_i = (l_xi, l_yi, l_zi), i = 1, …, n, be the directions on corresponding to each triangle shape, and we assume that l₁, …, l_n are a random sample from the Fisher–von Mises distribution. From the data we find that

Thus the resultant mean length is The mean directional vector is (0.146, 0.076, 0.472) and the resulting estimate of κ from the Fisher–von Mises distribution is , so we have a large value of the concentration here. When mapped to Bookstein coordinates the mean shape estimate is u^B = −0.101 and v^B = 0.632.

10.2.9 Simulation

There are several different methods available for simulation from the complex Bingham distribution, and Kent et al. (2004) provide several algorithms, which have rather different efficiencies depending on the amount of concentration. A detailed description of simulation methods for a variety of directional distributions is given by Kent et al. (2013). Also, Kume and Walker (2009) discuss simulation for the Fisher–Bingham distribution.

10.3.1 The density

Definition 10.2 The complex Watson distribution on the pre-shape sphere , denoted by , has p.d.f.

(10.9)

Properties of the distribution were discussed by Mardia and Dryden (1999). In shape analysis we are primarily interested in ξ ≥ 0, although the concentration parameter ξ can be negative in some circumstances. For ξ > 0, the modal pre-shape direction is μe^iθ, where 0 ≤ θ < 2π is an arbitrary rotation angle. The complex Watson distribution is an important special case of the complex Bingham distribution, when there are just two distinct eigenvalues in A (a single distinct largest eigenvalue and all other eigenvalues being equal). In this case all directions orthogonal to the modal axis have equal weight, and so this model implicitly assumes a spherical error distribution, as would be obtained from independent isotropic landmarks with equal variances. In this case the complex Bingham distribution can be re-parameterized so that A = −2κ(I_{k − 1} − μμ*), where μ is the modal vector on the pre-shape sphere. The density on the pre-shape sphere can be written as:

(10.10)

where ρ(H^Tz, H^Tμ) is the Procrustes distance between the shapes of H^Tz and H^Tμ. Note that H^Tz is the centred icon corresponding to the pre-shape z. The parameters of the distribution are the population mode pre-shape direction μ and the concentration parameter ξ = 2κ. We have chosen to re-parameterize the complex Watson distribution with the parameter κ = ξ/2 in order to make connections with the offset normal distributions below in Section 11.1.

We call this distribution the complex Watson distribution due to the similar form of the density to the Dimroth–Watson distribution on the real sphere, see Watson (1983, 1965), who called it the Scheiddegger–Watson distribution.

For triangles (k = 3) the complex Watson distribution is the same as the complex Bingham distribution, as there are at most two distinct eigenvalues in that case. An alternative expression for the density is

where d_F( · ) is the full Procrustes distance of Equation (4.10).

This distribution was briefly discussed by Dryden (1991), Kent (1994) and Prentice and Mardia (1995). Full details and further developments can be found in Mardia and Dryden (1999), on which subsequent details are based.

Result 10.6 The integrating constant for the complex Watson distribution is:

where

(10.11)

is the confluent hypergeometric function (e.g. see Abramowitz and Stegun 1970).

See Mardia and Dryden (1999) for a proof. Note that the integrating constant has a relatively simple form, and so inference is quite straightforward.

10.3.2 Inference

Let z₁, …, z_n be a random sample from a population modelled by a complex Watson distribution, with n ≥ k − 1. Set

to be the (k − 1) × (k − 1) complex sum of squares and products matrix. Suppose that the eigenvalues of S are positive with eigenvalues l_{k − 1} > l_{k − 2} ≥ … ≥ l₁ > 0 and let g_{k − 1}, …, g₁ denote the corresponding eigenvectors.

Result 10.7 Under the complex Watson distribution the MLEs of μ and κ are given by:

and the solution to

where Under high concentrations

(10.12)

Proof: The log-likelihood for the data reduces to:

Holding κ constant, it can be seen that provides the maximum, where α is an arbitrary rotation angle (0 ≤ α < 2π), and so

The MLE of κ is found by solving

Under high concentration

and so

Therefore

and so the high concentration approximate MLE for κ follows.

As for the complex Bingham distribution the dominant eigenvector can be regarded as the average axis of the data (up to a rotation), and it has the same shape as the sample full Procrustes mean of Equation (6.11).

Since the complex Watson distribution is a complex Bingham distribution with all but the largest eigenvalue equal, then under the complex Watson distribution we should expect l₁, …, l_{k − 2} to be approximately equal. In other words, the shape variability for the complex Watson distribution is isotropic and all principal components of shape variability have equal weightings. This observation provides us with a model checking procedure.

10.3.3 Large concentrations

We now consider the case of large κ.

Result 10.8 If κ is large, then the distribution tends to a complex multivariate normal distribution in the tangent space with mean 0 and covariance matrix , which is a generalized inverse of 2κ(I_{k − 1} − μμ*).

Proof: We show that the following expression is true:

(10.13)

where z*z = μ*μ = 1. Furthermore, A⁻ = A where A = I − μμ*, and A⁻ is a generalized inverse. First we prove the generalized inverse, that is AA⁻A = A. It is easily verified that A is idempotent and therefore A³ = A. Hence, a generalized inverse of A is also A. To prove Equation (10.13) expand the left-hand side:

Since z*z = 1, we have

Furthermore, (I − μμ*)μ = μ − μ(μ*μ) = μ − μ = 0. Hence Equation (10.13) is proved. Thus, the exponential term in the complex Watson density reduces to:

where v are the tangent plane coordinates, and z is close to μ. Hence,

for large concentrations.

For large concentrations the distribution is also very similar to the offset isotropic normal distribution described in Section 11.1.2. It follows from the normal approximation that for large κ

is approximately distributed as χ²_{2k − 4}, since the complex rank of (I − μμ*) is k − 2. Furthermore, we have

where ρ is the Riemannian distance between the shapes of H^Tz and H^Tμ, so that

(10.14)

This result agrees with the approximate distribution for the squared full Procrustes distance under isotropic normal errors given in Equation (9.6), with

where S(H^Tμ) is the centroid size of the icon H^Tμ corresponding to μ.

Important point: The central role that the Fisher–von Mises distribution plays in directional data analysis is played by the complex Watson distribution for 2D shape analysis. For triangles we have seen that the two distributions are equivalent. For more than three points there are also similarities. Both distributions are very tractable, and each can be viewed as a conditional normal distribution.

10.9.1 Rotationally symmetric size-and-shape family

We could specify size-and-shape distributions directly in size-and-shape space SΣ^k_m. A suitable candidate could be that based on Riemannian distance, with density proportional to:

where d_S( · ) was given in Equation (5.5) and ϕ( · ) ≥ 0 is a suitable penalty function. In the case ϕ(x) = x² the MLE of size-and-shape will be the same as those using least squares partial Procrustes analysis, provided the integrating constant does not depend on μ, that is for the planar case. As for the shape case in Section 10.6 we can choose different functions for estimators with different properties. Further size-and-shape distributions based on normally distributed data are discussed in Section 11.5.

10.9.2 Central complex Gaussian distribution

For 2D data we could exploit the complex representation and, in an analogous manner to using the complex Bingham distribution in pre-shape space for shape analysis, we could specify a distribution for the Helmertized landmarks which is invariant under rotations. A suitable candidate is the central complex Gaussian distribution,

where Σ is Hermitian. The distribution of z_H is identical to that of e^iθz_H for any arbitrary rotation θ. The mode of the distribution is given by the dominant eigenvector of Σ, hence the MLE of size-and-shape under this model is the dominant eigenvector of

given a random sample of Helmertized observations z_H1, …, z_Hn.