Fitting multivariate skew-normal distributions

Usage

msn_fit(X, y, freq, plotit, traceout, iter_max, x_tol )

Arguments

y a matrix or a vector. In y is a matrix, its rows refer to observations, and its columns to components of the multivariate distribution. In y is a vector, it is converted to a one-column matrix, and a scalar skew-normal distribution is fitted.
X a matrix of covariate values. If missing, a one-column matrix of 1's is created; otherwise, it must have the same number of rows of y.
freq a vector of weights. If missing, a one-column matrix of 1's is created; otherwise it must have the same number of rows of y.
plotit logical value which controls the graphical output (default=1); see below for description.
traceout logical value which controls printing of the algorithm convergence. If traceout=1, details are printed. Default value is 0.
iter_max maximum number of iterations in the maximisation routine (default is 150).
x_tol tolerance (default is 1e-8).
start starting values for the optimisation routine (see msn_mle).

Description

Fits a multivariate skew-normal (MSN) distribution to data, or fits a linear regression model with multivariate skew-normal errors, using maximum likelihood estimation. The outcome is then displayed in graphical form.

Details

For computing the maximum likelihood estimates, msn_fit invokes msn_mle which does the actual computational work; then, msn_fit displays the results in graphical form. The documentation of msn_mle gives details of the numerical procedure for maximum likelihood estimation.

Although the function accepts a vector y as input, the use of sn_mle is recommended in the scalar case.

Value

a list containing the following components:

dp a list containing the direct parameters beta, Omega, alpha. Here, beta is a matrix of regression coefficients with size(beta)=[size(X,1),size(y,2)], Omega is a covariance matrix of order size(y,2), alpha is a vector of shape parameters of length size(y,2).
logL log-likelihood evaluated at dp.
se a list containing the components beta, alpha, info. Here, beta and alpha are the standard errors for the corresponding point estimates; info is the observed information matrix for the working parameter, as explained below.
options the list returned by the optimizer routine; see the documentation of foptions for explanation of its components.
test_normality a list of with elements test and p_value, which are the value of the likelihood ratio test statistic for normality (i.e. test that all components of the shape parameter are 0), and the corresponging p-value.

Side Effects

Graphical output is produced if (missing(freq))=1. Three plots are produced, and the programs pauses between each two of them, waiting for any key to be pressed.

The first plot uses the variable y if X is missing, otherwise it uses the residuals from the regression. The form of this plot depends on the value of k=size(y,2); if k=1, an histogram is plotted with the fitted distribution suerimposed. If k>1, a matrix of scatterplots is produced, with superimposed the corresponging bivariate densities of the fitted distribution.

The second plot has two panels, each representing a QQ-plot of Mahalanobis distances. The first of these refers to the fitting of a multivariate normal distribution, a standard statistical procedure; the second panel gives the corresponding QQ-plot of suitable Mahalanobis distances for the multivariate skew-normal fit.

The third plot is similar to the previous one, except that PP-plots are produced.

BACKGROUND

The multivariate skew-normal distribution is discussed by Azzalini and Dalla Valle (1996); the (Omega,alpha) parametrization adopted here is the one of Azzalini and Capitanio (1998).

References

Azzalini, A. and Dalla Valle, A. (1996). The multivariate skew-normal distribution. Biometrika 83, 715-726.

Azzalini, A. and Capitanio, A. (1999). Statistical applications of the multivariate skew-normal distribution. J.Roy.Statist.Soc. B 61, part 3.

See Also

msn_mle, sn_mle

Examples

# a simple case
msn_fit(NaN,[lbm,bmi,ssf])  #no matrix X
#
# a regression case
a = msn_fit([ones(length(lbm),1),lbm], bmi, NaN, NaN, NaN, NaN, 1e-6)
# uses x_tol = 1e-6 and default values for the other input parameters
# refine the previous outcome
a1 = msn_fit([ones(length(lbm),1),lbm], bmi,  NaN, NaN, NaN, NaN, 1e-9,a_dp)


[Package Contents]