# Underlying Theory of the Test

## Sources

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.

## Transformations and the Characterization

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 support region for S

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~.

pos.def.plot <- function(npoints, v12){
nv12 <- length(v12)
x <- rep(1:npoints,each=npoints) /npoints
y <- rep(1:npoints,times=npoints) /npoints
v12 <- rep(v12,each=npoints*npoints)
nameV12 <- paste("|V12| =",v12)
calcpd <- x*y - v12*v12
posdef <- calcpd > 0
xdf <- data.frame(x,y,v12,calcpd,posdef,nameV12)
xdf[xdf$posdef,] } tempdf2 <- pos.def.plot(npoints=40,v12=c(.15,.30,.45,.60,.75,.9)) horlabel="V1" vertlabel="V2" legendname=NULL cap=NULL mtitle="Support for MVNchar Test" stitle="Positive Definite Region for p=2" XVAR <- tempdf2[,1] YVAR <- tempdf2[,2] COV1 <- tempdf2[,5] dfplot <- data.frame(XVAR,YVAR,COV1) FACET <- tempdf2[,6] dfplot <- data.frame(dfplot,FACET) out <- ggplot2::ggplot(data=dfplot,ggplot2::aes(XVAR,YVAR,COV1)) + ggplot2::geom_point() out$labels$shape <- legendname names(FACET)<-tempdf2$namesV12
out <- out + ggplot2::facet_wrap(~FACET)
out <- out + ggplot2::ggtitle(mtitle,subtitle=stitle) + ggplot2::xlab(horlabel) + ggplot2::ylab(vertlabel) +                                  ggplot2::labs(caption=cap,legend=legendname)
out <- out + ggplot2::labs(color=legendname)
print(out)


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.

# Implementing the Test

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.

## Relating Sample Size to the Size of Mini-cubes

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
r
a <- c(2,5,10,15,20,256,4000,31997,107989,255974,427,6666,53328,179982,426624)
a <- c(a,2,5,6,7,9,3648, 429504,1402224,3537024,15917712,6080,715840,2337040,5395040,26529520)
a <- c(a,2,3,4,5,6,20640,1224600,20132640,187000000,1190000000,34400,2041000,33554400,312000000,1980000000)
a <- matrix(a,nrow=5,byrow=FALSE)
b <- matrix("   ",nrow=5,ncol=1)
c <- matrix("         ",nrow=5,ncol=1)
a <- cbind(a,b,c)
a <- as.data.frame(a, row.names=NULL, optional=TRUE, col.names=NULL)
a <- a[,c(11,1:3,10,4:6,10,7:9)]
a



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.

# Example Databases Included in MVNtestchar Package

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}

### Theorem 1.

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~(μ,Σ).

### Proof.

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~.

### Theorem 2.

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 Σ .

### Proof.

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.

## Try the MVNtestchar package in your browser

Any scripts or data that you put into this service are public.

MVNtestchar documentation built on July 26, 2020, 1:07 a.m.