FAQ/Rmardia - CBU statistics Wiki

Upload page content

You can upload content for the page named below. If you change the page name, you can also upload content for another page. If the page name is empty, we derive the page name from the file name.

File to load page content from
Page name
Comment
Finzd thee wrang lelters ino eacuh wosrd

location: FAQ / Rmardia

MATLAB code for computing Mardia's multivariate kurtosis coefficient in Multivariate Normality testing Normality testing

David Graham (2006) has written MATLAB code based upon the work of A. Trujillo-Oriz and R. Hernandez-Walls (2003) from the Universidad Autonoma de Baja California who compute Mardia's multivariate kurtosis coefficient which, including an example, is located here.

Reference

Trujillo-Ortiz, A. and R. Hernandez-Walls. (2003). Mskekur: Mardia's multivariate skewness and kurtosis coefficients and its hypotheses testing. A MATLAB file. [WWW document]. URL http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?

This code is also reproduced below:

function [H stats] = mardiatest(X,alpha)
% MARDIATEST Mardia multivariate normality test using skewness & kurtosis.
%   H = MARDIATEST(X,ALPHA) performs Mardia's test of multivariate
%   normality to determine if the null hypothesis of multivariate 
%   normality is a reasonable assumption regarding the population
%   distributions of a random samples contained within the columns of X. X
%   must be a N-values by M-samples array. The desired significance level,
%   ALPHA, is an optional scalar input (default = 0.05).
%
%   The function performs three tests: of the multivariate skewness; the
%   multivariate skewness corrected for small samples; and the multivariate
%   kurtosis. H is a 1-by-3 array indicating the results of the hypothesis
%   tests:
%     H(i) = 0 => Do not reject the null hypothesis at significance level
%     ALPHA.
%     H(i) = 1 => Reject the null hypothesis at significance level ALPHA.
%
%   [H STATS] = MARDIATEST(...) also returns a structure array with the
%   following fields:
%     stats.Hs  - logical scalar indicating whether to reject the null
%                 hypothesis that the multivariate skewness is consistent
%                 with a multivariate normal distribution.  
%     stats.Ps  - the P value for the multivariate skewness.
%     stats.Ms  - the Mardia test statistic for the multivariate skewness.
%     stats.CVs - the critical value for the multivariate skewness.
%     stats.Hsc - logical scalar indicating whether to reject the null
%                 hypothesis that the multivariate skewness (corrected for
%                 small samples) is consistent with a multivariate normal
%                 distribution.
%     stats.Psc - the P value for the multivariate skewness (corrected for
%                 small samples).
%     stats.Msc - the Mardia test statistic for the multivariate skewness
%                 (corrected for small samples).
%     stats.Hk  - logical scalar indicating whether to reject the null
%                 hypothesis that the multivariate kurtosis is consistent
%                 with a multivariate normal distribution.  
%     stats.Pk  - the P value for the multivariate kurtosis.
%     stats.Mk  - the Mardia test statistic for the multivariate kurtosis.
%     stats.CVk - the critical value for the multivariate kurtosis.
%
%   Notes:
%   For multivariate data, tests of normality of the individual variables
%   are not sufficent to determine multivariate normality. Even if every
%   variable in a set is normally distributed, it is still possible that
%   the combined distribution is not multivariate normal. Mardia's test is
%   one means testing for multivariate normality. The test is based on
%   independent tests of multivariate skewness and kurtosis. Data can be
%   assumed to conform to multivariate normality only if the null
%   hypothesis is not rejected for both tests.
%
%   Example:
%   >> X = [2.4  2.1  2.4; ...
%           3.5  1.8  3.9; ...
%           6.7  3.6  5.9; ...
%           5.3  3.3  6.1; ...
%           5.2  4.1  6.4; ...
%           3.2  2.7  4.0; ...
%           4.5  4.9  5.7; ...
%           3.9  4.7  4.7; ...
%           4.0  3.6  2.9; ...
%           5.7  5.5  6.2; ...
%           2.4  2.9  3.2; ...
%           2.7  2.6  4.1];
% 
%   >> [H stats] = mardiatest(X, 0.05)
% 
%   H =
%        0     0     1
% 
%   stats = 
%        Hs: 0
%        Ps: 0.9660
%        Ms: 3.5319
%       CVs: 18.3070
%       Hsc: 0
%       Psc: 0.8918
%       Msc: 4.9908
%        Hk: 1
%        Pk: 0.0452
%        Mk: -1.6936
%       CVk: 1.6449
%   
%   The magnitude of the Mardia kurtosis test statistic (Mk) is greater
%   than critical value (CVk) for a 5% level test, so the hypothesis of
%   multivariate normality is rejected. 
%
%   See also JBTEST, KSTEST, KSTEST2, LILLIETEST.
%   
%   This version:
%   Copyright 2006 David Graham, Loughborough University
%   (D.J.Graham@lboro.ac.uk)

% This function is based on MSKEKUR by A. Trujillo-Ortiz and R.
% Hernandez-Walls. Modifications include:
%   - additional checking of inputs
%   - modification of help text for consistency with Matlab conventions
%   - the output of the test results to variables
%   - removal of the code sending test results as a string to the command
%     window
%   - removed call to 'eval'.
%
% Statements from the original code:
%
% %  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
% %         Facultad de Ciencias Marinas        
% %         Universidad Autonoma de Baja California 
% %         Apdo. Postal 453  
% %         Ensenada, Baja California
% %         Mexico  
% %         atrujo@uabc.mx
% %
% %  May 22, 2003.
% %
% %  To cite this file, this would be an appropriate format:
% %  Trujillo-Ortiz, A. and R. Hernandez-Walls. (2003). Mskekur: Mardia's multivariate skewness
% %    and kurtosis coefficients and its hypotheses testing. A MATLAB file. [WWW document]. URL 
% %    http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=3519&objectType=FILE
% %
% %  References:
% % 
% %  Mardia, K. V. (1970), Measures of multivariate skewnees and kurtosis with
% %         applications. Biometrika, 57(3):519-530.
% %  Mardia, K. V. (1974), Applications of some measures of multivariate skewness
% %         and kurtosis for testing normality and robustness studies. Sankhyâ A,
% %         36:115-128
% %  Stevens, J. (1992), Applied Multivariate Statistics for Social Sciences. 2nd. ed.
% %         New-Jersey:Lawrance Erlbaum Associates Publishers. pp. 247-248.


% Check number of input arguements
if (nargin < 1) || (nargin > 2)
    error('Requires one or two input arguments.')
end

% Define default ALPHA if only one input is provided
if nargin == 1, 
    alpha = 0.05;
end

% Check for validity of ALPHA
if ~isscalar(alpha) || alpha>0.5 || alpha <0
    error('Input ALPHA must be a scalar between 0 and 0.5.')
end

[n,p] = size(X);

% Check for validity of X
if p < 2 || ~isnumeric(X)
    error('Input X must be a numeric array with at least 2 columns.')
end


difT = [];
for     j = 1:p
   difT = [difT,(X(:,j) - mean(X(:,j)))];
end;

S = cov(X);                     % Variance-covariance matrix
D = difT * inv(S) * difT';      % Mahalanobis' distances matrix
b1p = (sum(sum(D.^3))) / n^2;   % Multivariate skewness coefficient
b2p = trace(D.^2) / n;          % Multivariate kurtosis coefficient

k = ((p+1)*(n+1)*(n+3)) / ...
    (n*(((n+1)*(p+1))-6));      % Small sample correction
v = (p*(p+1)*(p+2)) / 6;        % Degrees of freedom
g1c = (n*b1p*k) / 6;            % Skewness test statistic corrected for small sample (approximates to a chi-square distribution)
g1 = (n*b1p) / 6;               % Skewness test statistic (approximates to a chi-square distribution)
P1 = 1 - chi2cdf(g1,v);         % Significance value of skewness
P1c = 1 - chi2cdf(g1c,v);       % Significance value of skewness corrected for small sample

g2 = (b2p-(p*(p+2))) / ...
    (sqrt((8*p*(p+2))/n));      % Kurtosis test statistic (approximates to a unit-normal distribution)
P2 = 1-normcdf(abs(g2));        % Significance value of kurtosis


% Prepare outputs

stats.Hs  = (P1 < alpha);
stats.Ps  = P1;
stats.Ms  = g1;
stats.CVs = chi2inv(1-alpha,v);
stats.Hsc = (P1c < alpha);
stats.Psc = P1c;
stats.Msc = g1c;
stats.Hk  = (P2 < alpha);
stats.Pk  = P2;
stats.Mk  = g2;
stats.CVk = norminv(1-alpha,0,1);

H = [stats.Hs stats.Hsc stats.Hk];