MVMMT
Computes Mardia’s multivariate measures of skewness and kurtosis and test for multivariate normality.
Required Arguments
X — NOBS by NVAR+ m matrix containing the data. (Input)
m is 0, 1, or 2 depending upon whether any columns in X contain frequencies or weights.
IND — Vector of length NVAR containing the column numbers in X for which statistics are desired. (Input)
STAT — Vector of length 13 containing the output statistics. (Output)
If a statistic is not computed, the corresponding element of STAT is set to not a number (NaN).
STAT(1) = estimated skewness.
STAT(2) = expected skewness assuming a multivariate normal distribution.
STAT(3) = asymptotic chi‑squared statistic assuming a multivariate normal distribution.
STAT(4) = probability of a greater chi‑squared.
STAT(5) = Mardia and Foster’s standard normal score for skewness.
STAT(6) = estimated kurtosis.
STAT(7) = expected kurtosis assuming a multivariate normal distribution.
STAT(8) = asymptotic standard error of the estimated kurtosis.
STAT(9) = standard normal score obtained from STAT(6) through STAT(8).
STAT(10) = p‑value corresponding to STAT(9).
STAT(11) = Mardia and Foster’s standard normal score for kurtosis.
STAT(12) = Mardia’s SW statistic based upon STAT(5) and STAT(11).
STAT(13) = p‑value for STAT(12).
STAT(12) and STAT(13) are only computed when ICMPUT = 0.
Optional Arguments
NOBS — Number of rows of data in X. (Input)
Default: NOBS = size (X,1).
NVAR — Dimensionality of the multivariate space for which the skewness and kurtosis are to be computed. (Input)
Default: NVAR = size (IND,1).
NCOL — Number of columns in matrix X. (Input)
Default: NCOL = size (X,2).
LDX — Leading dimension of X exactly as specified in the dimension statement in the calling program. (Input)
Default: LDX = size (X,1).
IFRQ — Frequency option. (Input)
IFRQ = 0 means that all frequencies are 1.0. Positive IFRQ indicates that column number IFRQ of X contains the frequencies. All frequencies should be integer values. The NINT (nearest integer) function is used to obtain integer frequencies if this is not the case.
Default: IFRQ = 0.
IWT — Weighting option. (Input)
IWT = 0 means that all weights are 1.0. Positive IWT means that column IWT of X contains the weights. Negative weights are not allowed.
Default: IWT = 0.
ICMPUT — Option parameter giving the statistics to compute. (Input)
Default: ICMPUT = 0.
ICMPUT | Output Statistics |
|---|
0 | Both skewness and kurtosis. |
1 | Kurtosis only. |
2 | Skewness only. |
NI — The sum of the frequencies of all observations used in the computations. (Output)
SWT — The sum of the weights times the frequencies for all observations used in the computations. (Output)
XMEAN — Vector of length NVAR containing the sample means. (Output)
R — NVAR by NVAR upper triangular matrix containing the Cholesky RTR factorization of the covariance matrix. (Output)
LDR — Leading dimension of R exactly as specified in the dimension statement in the calling program. (Input)
NRMISS — Number of rows of data in X containing any missing values (NaN, not a number). (Output)
Rows with missing values in the columns IND, IFRQ, and IWT are excluded from the analysis.
FORTRAN 90 Interface
Generic: CALL MVMMT (X, IND, STAT [, …])
Specific: The specific interface names are S_MVMMT and D_MVMMT.
FORTRAN 77 Interface
Single: CALL MVMMT (NOBS, NVAR, NCOL, X, LDX, IND, IFRQ, IWT, ICMPUT, NI, SWT, XMEAN, R, LDR, STAT, NRMISS)
Double: The double precision name is DMVMMT.
Description
Routine MVMMT computes Mardia’s (1970) measures b1,p and b2,p of multivariate skewness and kurtosis, respectfully, for p = NVAR. These measures are then used in computing tests for multivariate normality. Three test statistics, one based upon b1,p alone, one based upon b2,p alone, and an omnibus test statistic formed by combining normal scores obtained from b1,p and b2,p are computed. On the order of np3, operations are required in computing b1,p when the method of Isogai (1983) is used, where n = NOBS. On the order of np2, operations are required in computing b2,p.
Let
where
fi is the frequency of the i‑th observation, and wi is the weight for this observation. (Weights wi are defined such that xi is distributed according to a multivariate normal, N(μ, Σ/wi) distribution, where Σ is the covariance matrix.) Mardia’s multivariate skewness statistic is defined as:
while Mardia’s kurtosis is given as:
Both measures are invariant under the affine (matrix) transformation AX + D, and reduce to the univariate measures when p = NVAR = 1. Using formulas given in Mardia and Foster (1983), the approximate expected value, asymptotic standard error, and asymptotic p‑value for b2,p, and the approximate expected value, an asymptotic chi‑squared statistic, and p‑value for the b1,p statistic are computed. These statistics are all computed under the null hypothesis of a multivariate normal distribution. In addition, standard normal scores W1(b1,p) and W2(b2,p) (different from but similar to the asymptotic normal and chi‑squared statistics above) are computed. These scores are combined into an asymptotic chi‑squared statistic with two degrees of freedom:
This chi‑squared statistic may be used to test for multivariate normality. A p‑value for the chi‑squared statistic is also computed.
Comments
1. Workspace may be explicitly provided, if desired, by use of M2MMT/DM2MMT. The reference is:
CALL M2MMT (NOBS, NVAR, NCOL, X, LDX, IND, IFRQ, IWT, ICMPUT, NI, SWT, XMEAN, R, LDR, STAT, NRMISS, D, OB, CC)
The additional arguments are as follows:
D — Work vector of length NVAR.
OB — Work vector of length NVAR.
CC — Work vector of length m, where m = NVAR * NVAR if ICMPUT = 1 or
m = NVAR * NVAR * NVAR otherwise.
2. Informational errors
Type | Code | Description |
|---|
4 | 1 | At least one of the variables in X is linearly related to the other variables in X. |
4 | 2 | The sum of the frequencies must be greater than the maximum of 3 and the number of variables plus one. |
Example
In the following example, 150 observations from a 5 dimensional standard normal distribution are generated via routine
RNNOR (see
Chapter 18, “Random Number Generation”). The skewness and kurtosis statistics are then computed for these observations.
USE IMSL_LIBRARIES
INTEGER LDR, LDX, NCOL, NVAR, I
PARAMETER (NCOL=5, LDX=150, NVAR=NCOL, LDR=NVAR)
!
INTEGER IND(5), NI, NOUT, NRMISS
REAL R(LDR,NVAR), STAT(13), SWT, X(LDX,NCOL), XMEAN(NVAR)
!
DATA IND/1, 2, 3, 4, 5/
!
CALL RNSET (123457)
DO 10 I=1, NCOL
CALL RNNOR(X(:, I))
10 CONTINUE
!
CALL MVMMT (X, IND, STAT, NI=NI, SWT=SWT, XMEAN=XMEAN, R=R, &
NRMISS=NRMISS)
!
CALL UMACH (2, NOUT)
WRITE (NOUT,*) ' NI = ', NI, ' SWT = ', SWT, ' NRMISS = ', NRMISS
CALL WRRRN ('XMEAN', XMEAN, 1, NVAR, 1)
CALL WRRRN ('R', R)
CALL WRRRN ('STAT', STAT, 1, 13, 1)
!
END
Output
NI = 150 SWT = 150.0 NRMISS = 0
XMEAN
1 2 3 4 5
0.0355 0.0467 0.0599 0.0957 0.1007
R
1 2 3 4 5
1 1.033 -0.022 -0.037 0.055 -0.003
2 0.000 0.993 -0.119 -0.076 -0.056
3 0.000 0.000 0.997 -0.089 0.017
4 0.000 0.000 0.000 1.008 -0.040
5 0.000 0.000 0.000 0.000 1.027
STAT
1 2 3 4 5 6 7 8 9 10
1.52 1.36 38.71 0.31 0.42 34.21 34.54 1.27 -0.26 0.80
11 12 13
0.18 0.21 0.90
