This test of multivariate normality is based on the univariate work of Csörgö and Seshadri (1970, 1971) and on my doctoral dissertation (Fairweather 1973). As indicated in the titles of these articles, the test is based on a characterization; that is, it is based on characteristics of the normal distribution that are unique to it among all (nondegenerate) multivariate distributions. The MVNtestchar package implements this test.

Consider a set of p x 1 random vectors of full rank X_{i}, i = 1, . . ., 4(p+1). Let Y_{i} = X_{2i} – X_{2i-1}, i = 1, . . ., 2(p+1).

The Y_{i} have a distribution that is centered at 0. In fact, all of the odd moments of the Y_{i} are zero regardless of the underlying distribution of the X_{i}. Now, define

\[ W_{1} = \sum_{i = 1}^{p + 1}{Y_{i}{Y'}_{i}} \]

and

\[ W_{2} = \sum_{i = p + 2}^{2(p + 1)}{Y_{i}{Y'}_{i}} \],

where Y’_{i} is the transpose of Y_{i}. The W_{i} are then independently distributed, symmetric matrices of rank p.

Let T = W_{1} + W_{2} and let S = T^{-1/2} W_{1} T’^{-1/2} . S is a positive definite (symmetric) matrix of rank p regardless of the underlying distribution of the X_{i}.

It is shown in the Appendix that S is distributed uniformly on its support region **if and only if** the X_{i} are multivariate normal. It is this characteristic that underlies the test.

The p x p symmetric matrix S is equivalent to a set of p(p+1)/2 random variables V_{1}, V_{2}, … , V_{p}, V_{12}, V_{13}, .., V_{1p}, …, V_{p-1,p}. This is easily seen if we lay out the V_{i} and the V_{ij} in the matrix format, showing only the upper triangle:

\[ \begin{matrix} \begin{matrix} V_{1} & V_{12} \\ & V_{2} \\ \end{matrix} & \begin{matrix} V_{13} & \begin{matrix} \cdots & V_{1,p} \\ \end{matrix} \\ V_{23} & \begin{matrix} \cdots & V_{2,p} \\ \end{matrix} \\ \end{matrix} \\ \begin{matrix} & \\ & \\ \end{matrix} & \begin{matrix} \begin{matrix} V_{3} & \cdots \\ \end{matrix} & V_{3,p} \\ \begin{matrix} & \\ \end{matrix} & V_{p} \\ \end{matrix} \\ \end{matrix} \]

V_{1} through V_{p} are the diagonal elements of the matrix and V_{12} … V_{p-1,p} are the off-diagonal elements.

By construction, the diagonal elements of S all lie in the interval [0,1]. Because S is positive definite, the off-diagonal elements all lie in the interval [-1,1]. The support region of S is within the hyper-rectangle

\[ R_p = \lbrack 0,1 \rbrack ^p \lbrack -1,1 \rbrack ^m \]

where m = p(p-1)/2.

The support region is difficult to envision in higher dimensions. For p = 2, the matrix S has two diagonal elements and one off-diagonal element: V_{1}, V_{2}, and V_{12}. The support region of S is then that part of the 3-dimensional rectangle bounded by R_{2} = [0,1] x [0,1] x [-1,1] that satisfies V_{1} V_{2} – V_{12}^{2} > 0.

By stepping |V_{12}| up from 0 to 1.0, the following graph shows the support region (shaded) as a function of V_{1} and V_{2}.

The package contains four functions that produce rotatable 3-dimensional graphs depicting this hyper-rectangle and the positive definite region within it. The function *support.p2( )* shows the entire positive definite region. The functions *slice.v1( )* and *slice.v12( )* show slices through the region for fixed values of V_{1} and V_{12}, respectively. Finally, *maxv12( )* shows the maximum value of V_{12} as a function of V_{1} and V_{2}.

The latest values of the rotational parameters are output in list format upon exit from the graphics functions to facilitate return to that rotation, if desired. To capture these, assign the function to a new variable.

The function *testunknown(x, pvector, k)* implements the test of multivariate normality of the n x p matrix, x. The input parameters will be discussed more fully in the next paragraphs.

The matrix x is assumed to be a sample of n observations on the unknown p-variate distribution. Here, n = 4r(p+1) for r a positive integer. The transformations involve exact numbers of random variables, and *testunknown( )* will discard observations at random to ensure that this condition holds.

S_{1}, …, S_{r} are distributed uniformly on the positive definite subspace of the hyper-rectangle

\[ R_p~ = \lbrack 0,1 \rbrack ^p \lbrack -1,1 \rbrack ^m \]

if and only if x is a sample from a multivariate normal distribution.

*testunknown( )* performs a chisquare goodness of fit test by filling this hyper-rectangle with mini-cubes and counting and comparing the number of S_{i} in each mini-cube. As implied, the mini-cubes are hyper-cubical (equal length in every dimension).

The input variable pvector is essentially a check to ensure that the matrix is oriented properly; it should equal the value of p taken from x within the function. If this test fails, the function aborts.

The parameter k defines the number of cuts to be made on each edge of the hyper-rectangle. The mini-cubes will then be of length 1/k on each edge. There are p(p+1)/2 terms in this set.

It is possible to undertake various research projects with this test function, and an array with mobs=b layers is allowed in order to facilitate this possibility. We have in mind simulations with mobs repetitions. If x is an array with b layers, x should have dimension n x p x b.

*testunknown( )* divides the sample into r sets of 4(p+1) observations and performs the transformations described above on each of the r sets independently. This results in positive definite matrices S1, …, Sr distributed independently as described above regardless of the distribution of the X_{i}.

S_{1}, …, S_{r} are distributed uniformly on the positive definite subspace of the hyper-rectangle \(R_p = \lbrack 0,1 \rbrack ^p \lbrack - 1,1 \rbrack^m\) if and only if x is a sample from a multivariate normal distribution.

Any goodness of fit test, univariate or multivariate, must consider the relationship of the sample size to the number of “bins” into which the support is subdivided. The sample size is always finite and the samples are continuous, so that creating too many bins will always result in exactly one observation per occupied bin. With too many bins, there can be no distinction between null and alternative distributions.

In our case, the number of mini-cubes (“bins”) into which the hyper-cube is divided is N(k,p) = k^{p}(2k)^{m},where m = p(p-1)/2. For p = 2, m = 1 and N(k,p) = 2k^{3}. Similarly, N(k,3) = 8k^{6} and N(k,4) = 64k^{10} . Mini-cubes are entirely, partially, or not at all within the positive definite region of the hyper-rectangle. We can calculate analytically the fraction of the hyper-rectangle that is within the positive definite region only for p = 2. This fraction is 4/9. For p > 2, the calculation appears to be intractable.

The number of mini-cubes clearly grows rapidly with the dimension of the sample. However, the support of the S_{i} is only a subset of the hyper-rectangle, namely the positive definite region. The following table shows the number of mini-cubes in the hyper-rectangle, N(k,p), the ratio of the positive definite region to the overall volume of the hyper-rectangle, and the approximate number of mini-cubes in the positive definite region, for several values of k and p.

**Table 1. The number of mini-cubes, N(k,p) in the hyper-rectangle, as a function of the number of cuts, k and the dimensionality, p of the sample, and the number that represent positive definite matrices.**

```
-------------------------- -------------------------- --------------------------
p = 2 p=3 p=4
- ------ ----- ------- - ------- ----- ------- - -------- ----- -------
k N(k,p) ratio pos def k N(k,p) ratio pos def k N(k,p) ratio pos def
(%) (%) (%)
2 16 44.4 7 2 512 14.8 76 2 65536 0.6 344
5 250 44.4 111 5 125000 7.2 8948 3 3.7e6 0.5 20410
10 2000 44.4 889 6 373248 7.8 29213 4 6.7e7 0.5 335544
15 6750 44.4 3000 7 941192 7.8 73688 5 6.3e8 0.5 3.0e6
20 16000 44.4 7110 9 4.2e6 7.8 331619 6 3.9e9 0.5 2.0e7
```

N(k,p) is easily calculated in each case. For p=2, the calculated ratio was multiplied by the number of mini-cubes to get the approximate number of positive definite mini-cubes. For p > 2, rows 1 through 4 of Table 1 were calculated as follows: The hyper-rectangle was filled with mini-cubes as described above. A mini-cube was defined to be within the positive definite region if a point very near the center of the mini-cube represented a positive definite matrix. The last row of the table is an extrapolation obtained by applying the asymptotic ratio to the calculated value of N(k,p).

For each value of p the ratio of mini-cubes in the positive definite region to the overall number in the hyper-rectangle is fairly constant. For p > 2, this ratio is a very small part of the overall volume of the hyper-rectangle. Nevertheless, Table 1 shows that a very large number of “bins” in the support region will result if k is set too high.

In performing the characterization transformations, the number of vector samples is substantially reduced to form the positive definite matrices that are tested for uniformity of distribution. Table 1 refers to the number of bins into which the matrices S_{i} will fall. 4(p+1) vectors X_{i} will result in a single matrix S_{i}. This multiplier is 12 for p = 2, is 16 for p = 3, and is 20 for p = 4. Assuming that the expected number of S_{i} in each bin should be about E, Table 2 gives the number of X_{i} that should be in the sample for each value of k.

**Table 2. Relationship of sample size n to number of cuts k, as a function of the expected number E of positive definite matrices Si per mini-cube.**

```
--------------- -------------------- ---------------------
p = 2 p = 3 p = 4
- ---- ---- - ----- ---- - ---- ----
k E=3 E=5 k E=3 E=5 k E=3 E=5
```
##
## 1 2 256 427 2 3648 6080 2 20640 34400
## 2 5 4000 6666 5 429504 715840 3 1224600 2041000
## 3 10 31997 53328 6 1402224 2337040 4 20132640 33554400
## 4 15 107989 179982 7 3537024 5395040 5 1.87e+08 3.12e+08
## 5 20 255974 426624 9 15917712 26529520 6 1.19e+09 1.98e+09
```
```

From Table 2, we conclude that if we wish to test a trivariate sample for normality, and we have about 720,000 observations, we should set k to be no more than 5 to ensure that there are about E=5 observations per mini-cube.

In order to gain experience with the functions, the package contains four sample databases:

unknown.Np2

unknown.Np4

unknown.Bp2

unknown.Bp4

Here, bivariate vectors are symbolized by 2 and quadri-variate vectors are symbolized by 4. N symbolizes normal random vectors and B symbolizes modified Bernoulli random vectors. True Bernoulli random variables cause the test program to crash because of colinearity, so a normal variable with extremely small variance was added to make the Bernoulli vectors continuous random variables. unknown.Bp2 is a matrix; the others are arrays with a single layer.

*References*

Anderson, TW. (1958), An Introduction to Multivariate Statistical Analysis, John Wiley, New York.

Cramér, H (1962). Random Variables and Probability Distributions, Cambridge University Press, London.

Csörgö, M and Seshadri, V (1970). On the problem of replacing composite hypotheses by equivalent simple ones, *Rev. Int. Statist. Instit.*, **38**, 351-368

Csörgö,M and Seshadri,V (1971). Characterizing the Gaussian and exponential laws by mappings onto the unit interval, *Z. Wahrscheinlickhkeitstheorie verw. Geb.*, **18**, 333-339.

Deemer,WL and Olkin,I (1951). The Jacobians of certain matrix transformations useful in multivariate analysis, *Biometrika*, **58**, 345 367.

Fairweather, WR (1973). A test for multivariate normality based on a characterization. Dissertation submitted in partial fulfillment of the requirements for the Doctor of Philosophy, University of Washington, Seattle WA.

**APPENDIX**

The two theorems proven in this section characterize the multivariate normal probability distribution in terms of the distribution of certain positive definite matrix statistics. We will employ the following notation, definitions and conventions:

Capital letters will denote vectors or matrices, which should be clear by context or specific wording. A’ will denote the transpose of A. The determinant of a square matrix A will be denoted by |A| and its trace by tr A.

As usual, “independent and identially distributed” will be abbreviated as iid; “positive definite” as pd ; “is distributed as” by the symbol ~ ; and “is proportional to” by the symbol \(\propto\).

We write A < B if A and B are symmetric matrices and B – A is positive definite.

Define the functions

\[ \Gamma\left( z \right) = \int_{0}^{\infty}{t^{z - 1}\varepsilon^{- t}\text{dt}} \]

\[ \Gamma_{p}\left( \lambda \right) = \pi^{p(p - 1)/4}\prod_{j = 1}^{p}{\Gamma(\lambda - \frac{1}{2}}(j - 1)) \]

\[ \beta_{k,p}^{*}\left( a_{1},\ \ldots,\ a_{k} \right) = \frac{\left\{ \prod_{j = 1}^{k}{\Gamma_{p}\left( a_{j} \right)} \right\}}{\Gamma_{p}\left( \sum_{j = 1}^{k}a_{j} \right)} \]

\[ \beta_{k,p} = \beta_{k,p}^{*}(\frac{1}{2}\left( p + 1 \right),\ \ldots,\frac{1}{2}(p + 1)) \]

for \(\lambda\) and the a_{j} > (p-1)/2.

For S p.d., define

\[ c_{p}^{*}\left( n,\Sigma \right) = 1/\left\{ 2^{np/2}\Gamma_{p}\left( \frac{n}{2} \right){|\Sigma|}^{n/2} \right\} \]

\[ c_{p}\left( \Sigma \right) = c_{p}^{*}(p + 1,\Sigma). \]

The p-variate Normal distribution function with mean μ and covariance matrix Σ will be denoted by N_{p}(μ,Σ). The Wishart distribution with parameters Σ, n, and p will be denoted by W(Σ,n,p). If n ≥ p, a matrix variable S ~ W(Σ,n,p) has the density

\[ c_{p}^{*}\left( n,\Sigma \right){|S|}^{(n - p - 1)/2}exp( - \frac{1}{2}\text{trS}\Sigma^{- 1}), 0 < S. \]

We say that T is a matrix square root of U if TT’ = U. No further specification of T will be needed here.

The distance d(A,B) between symmetric p x p matrices A = ((a_{ij})) and B = ((b_{ij})) is defined to be

\[ d\left( A,B \right) = \left\lbrack \sum_{i = 1}^{p}{\sum_{j = 1}^{p}\left( a_{\text{ij}} - b_{\text{ij}} \right)^{2}} \right\rbrack^{1/2} \]

A function f of a symmetric p x p matrix is continuous at A if f is continuous at A in the metric d; f is continuous if it is continuous at every A. The function f is right continuous at A if for every ε > 0, all B > A for which d(A,B) < δ(ε) satisfy |f(A) – f(B)| < ε; f is right continuous if it is right continuous at every A.

Define the set of functions

V_{p} = {f: f defined on symmetric p x p matrices, f continuous at each A > 0 and right continuous at 0}.

G_{k,p} = {S_{i}: S_{i} is a p x p real, symmetric matrix, 0 < S_{1} < S_{2} < … < S_{k-1} < I}

Let Z_{i} (1 ≤ i ≤ 2mk) be iid p-variate random vectors with mean μ and pd covariance matrix Σ . Define

\[X_{i} = (Z_{2i-1} - Z_{2i})/ \sqrt{2}\] for (1 ≤ i ≤ mk) and

\[W_{j} = \sum_{i = \left( j - 1 \right)m + 1}^{\text{jm}}{X_{i}{X'}_{i}}\] for (1 ≤ j ≤ k).

Then the X_{i} are iid N_{p}(0,Σ) and the W_{j} are iid W(Σ,m,p) if and only if the Z_{i} ~ N_{p}(μ,Σ).

The X_{i} are iid by construction; similarly, for the W_{j}. The stated moments of the X_{i} also follow regardless of the distribution of the Z_{i}.

Sufficiency follows from the usual development of the Wishart distribution as given, for example, in Anderson (1958, Ch.7). We note that W_{j} does not necesarily have a density (i.e., may be degenerate) because the parameter m is not required here to be at least equal to p. The characteristic function of W_{j} does exist, however, in all cases.

To show that the condition is also necessary define B to be a non-singular matrix such that for θ real symmetric B’θB = D diagonal and also B’S^{-1}B = I. (Such a matrix exists; see Anderson (1958, p.34).) Then B^{-1}SB’^{-1} = I.

Let Y = B^{-1}X_{1}. The characteristic function of W_{1} for θ real symmetric is

\[ \psi_{W_{1}}\left( \theta \right) = \psi_{\sum_{i = 1}^{m}{X_{i}{X'}_{i}}}\left( \theta \right) = \left\lbrack \psi_{X_{i}{X'}_{i}}\left( \theta \right) \right\rbrack^{m} \]

because the X_{i} are iid. Now,

\[ \psi_{X_{i}{X'}_{i}}\left( \theta \right) = \psi_{BYY'B'}\left( \theta \right) \]

\[ = E\ exp\left( i\ tr\ BYY'B'\theta \right) \]

\[ = E\ exp\left( i\ Y'B'\theta BY \right) \]

\[ = E\ exp\left( i\ Y'DY \right) \]

\[ = E\ exp\left( \text{i}\sum_{j = 1}^{p}{d_{\text{jj}}Y_{j}^{2}} \right) \]

\[ = \psi_{S}\left( d \right) \]

This last term is the characteristic function of the vector S’ = (Y_{1}^{2}, …, Y_{p}^{2}), where d’ = (d_{11},…, d_{pp}) because d may be any vector with real elements by choosing θ appropriately. It follows that

\[ \psi_{W_{1}}\left( \theta \right) = \left\lbrack \psi_{S}(d) \right\rbrack^{m} \]

By assumption W_{1} ~ W(Σ,m,p); then for θ real symmetric (Anderson (1958, p.160)),

\[ \psi_{W_{1}}\left( \theta \right) = \frac{\left| \Sigma^{- 1} \right|^{m/2}}{\left| \Sigma^{- 1} - 2i\theta \right|^{m/2}} \]

\[ = \frac{\left( \left| B' \right|\left| \Sigma^{- 1} \right|\left| B \right| \right)^{m/2}}{\left( \left| B' \right|\left| \Sigma^{- 1} - 2i\theta \right|\left| B \right| \right)^{m/2}} \]

\[ = \left| I - 2iD \right|^{- m/2} \]

\[ = \left( \prod_{j = 1}^{p}\left( 1 - 2id_{\text{jj}} \right)^{- 1/2} \right)^{m} \]

This shows the components Y_{j}^{2} of S to be iid chi-square with 1 df. Write the j-th row of B^{-1} as b_{j}^{-1};

then Y_{j} = b_{j}^{-1}X_{1}. Because X_{1} has marginal symmetry about 0, Y_{j} is symmetrically distributed about 0 (1 ≤ j ≤ p). By Lemma 3 of Csörgö and Seshadri (1971), Y_{j} ~ N_{1}(0,1). Thus, Y ~ N_{p}(0,I) and X_{1} = BY ~ N_{p}(0,BB’), or X_{1} ~ N_{p}(0,S). The result follows by applying Cramér’s Theorem (Cramér, 1962, p. 112-113)) and the identity of distribution of the Z_{i}.

Let W_{i} (1 ≤ i ≤ k) be iid p x p pd random matrices with density f ε V_{p}. Define \[
U = \ \sum_{i = 1}^{p}W_{i}
\]

and let T be any matrix square root of U. Define U_{i} = T^{-1}W_{i}T’^{-1} and

\[ S_i = \sum_{j = 1}^{p}U_{j} \]

for (1 ≤ i ≤ k-1). Then (S_{1}, …, S_{k-1}) is uniform over G_{k,p} and independent of U if and only if the W_{i} ~ W(Σ,p+1,p) for some Σ .

Let W_{i} ~ W(Σ,p+1,p) (1 ≤ i ≤ k), so that

\[f\left( w_{i} \right) = c_{p}\left( \Sigma \right)\exp\left( -
\frac{1}{2}\text{tr}w_{i}\Sigma^{- 1} \right)\] for 0 < w_{i} .

Then because the Jacobian is unity,

\[ f_{w_{1},\ldots,w_{k - 1},u}\left( w_{1},\ldots,u \right) = \left\lbrack c_{p}(\Sigma) \right\rbrack^{k}\exp\left( - \frac{1}{2}\text{tru}\Sigma^{- 1} \right) \]

over \[\left\{ 0 < w_{i},\sum_{i = 1}^{k - 1}w_{i} < u \right\}\] .

The Jacobian of the transformation U_{i} = T^{-1}W_{i}T’^{-1} is shown by Deemer and Olkin (1951) to be |T^{-1}|^{p+1} . It follows that

\[ f_{U_{1},\ldots,U_{k - 1},U}\left( u_{1},\ldots,u \right) = \left( c_{p}\left( \Sigma \right) \right)^{k}\left| u \right|^{(k - 1)(p + 1)/2}\exp\left( - \frac{1}{2}\text{tru}\Sigma^{- 1} \right) \]

over {0 < u_{i}, \(\sum_{i = 1}^{k - 1}u_{i}\) < I, 0 < u}. The transformation from the U_{i} to the S_{i} has unit Jacobian and so

\[ f_{S_{1},\ldots,S_{k - 1},U}\left( s_{1},\ldots u \right) = \left( c_{p}\left( \Sigma \right) \right)^{k}\left| u \right|^{(k - 1)(p + 1)/2}\exp\left( - \frac{1}{2}\text{tru}\Sigma^{- 1} \right) \]

over {0 < s_{1} < … < s_{k-1} < I, 0 < u} . This density is seen to factor into the product of a constant and the W(S,k(p+1),p) density of U. We have thus demonstrated that U is independent of {S_{1}, …, S_{k-1}} and that

\[f_{S_{1},\ldots,S_{k - 1}}\left( s_{1},\ldots,s_{k - 1} \right) =
\beta_{k,p}^{- 1}\] for (s_{1}, …, s_{k-1}) ε G_{k,p}.

The converse is shown by developing a Cauchy functional equation for matrix arguments.

By assumption

\[ f_{S_{1},\ldots,S_{k - 1}|U} = f_{S_{1},\ldots,S_{k - 1}} = \beta_{k,p}^{- 1} \]

over G_{k,p} and for all U > 0. By construction

\[ f_{U_{1},\ldots,U_{k - 1},U}\left( u_{1},\ldots,u \right) = \left| u \right|^{- (k - 1)(p + 1)/2}f\left( u - \sum_{i = 1}^{k - 1}w_{i} \right)\prod_{i = 1}^{k - 1}{f\left( w_{i} \right)} \]

where u_{i} = s_{i} – s_{i-1} (s_{0}=0), w_{i} = tu_{i}t’ and tt’ = u. Then a density h of U satisfies

\[ h = \frac{f_{S_{1},\ldots,S_{k - 1},U}}{f_{S_{1},\ldots,S_{k - 1}|U}} \]

except possibly on a set of measure zero. Thus, for almost all (w_{1}, …, w_{k-1},u)

—-Equation 1—-:

\[ h\left( u \right) = \beta_{k,p}{|u|}^{- (k - 1)(p + 1)/2}f\left( u - \sum_{i = 1}^{k - 1}w_{i} \right)\prod_{i = 1}^{k - 1}{f(w_{i})}. \]

By the continuity assumption Equation 1 is satisfied at w_{1} = . . . = w_{k-1} = 0 for almost all u, which implies

\[ f\left( u - \sum_{i = 1}^{k - 1}w_{i} \right)\prod_{i = 1}^{k - 1}{f\left( w_{i} \right)} = f\left( u \right)f^{k - 1}\left( 0 \right)\text{ a.e.} \]

Therefore, f(0) > 0, Equivalently,

—-Equation 2—-:

\[ \prod_{i = 1}^{k}{g\left( w_{i} \right)} = g\left( \sum_{i = 1}^{k}w_{i} \right)\text{a.e.} \]

where \[w_{k} = u - \sum_{i = 1}^{k - 1}w_{i}\] and g = f^{-1}(0)f.

Equation 2 is a Cauchy functional equation for matrix arguments. In the present context g is a function of p(p+1)/2 (mathematically) independent variables, so that without loss of generality we may restrict consideration to upper triangular arguments. In particular, if W_{r} (1 ≤ r ≤ k) have all but the (i,j)-th element zero, Equation 2 is the usual Cauchy functional equation whose nontrivial solution is

—-Equation 3—-:

\[ g\left( W_{r} \right) = exp\left( - \frac{1}{2}a_{\text{ji}}w_{\text{ij}}^{\left( r \right)} \right)\text{} \]

With obvious notation, where a_{ji} is arbitrary.

Moreover, every upper triangular matrix W may be written as the sum of p(p+1)/2 matrices each of which has at most one non-zero element. Applying Equation 3 once again to this representation, we conclude that the most general solution to Equation 2 is

\[ g\left( W \right) = \prod_{i = 1}^{p}{\prod_{j = 1}^{p}{\exp\left( - \frac{1}{2}a_{\text{ji}}w_{\text{ij}} \right)}} \]

which allows g to be possibly non-trivial in each of its p(p+1)/2 arguments. In terms of symmetric matrices W,

\[ g\left( W \right) = exp\left( - \frac{1}{2}a_{\text{ji}}w_{\text{ij}} \right) \]

\[ = exp(- \frac{1}{2} tr WA) \]

where A is a real, symmetric matrix.

In order that f be a density proportional to g over the range 0 < W, one must have

\[ \sum_{i,j}^{}{w_{\text{ij}}a_{\text{ji}} > 0.} \]

(Otherwise, there exists a set of W’s with infinite Lebesgue measure in p(p+1)/2-space with densities greater than unity.) More specifically, by choosing w_{ij}’s appropriately, it can be shown that A must be pd.

We have just shown that W_{1} ~ W(A^{-1},p+1,p), which implies EW = (p+1)A^{-1}, and therefore A = Σ^{-1}, completing the proof.