Simon Preston

Analysis of spherical data with ESAG

This page contains a brief tutorial on how to use code (written jointly with Phillip Paine) which accompanies the paper:

An elliptically symmetric angular Gaussian distribution
Paine, P.J., Preston, S.P., Tsagris, M., Wood, A.T.A., Statistics and Computing (2017). doi:10.1007/s11222-017-9756-4 pdf d/l link

The code and data set used on this page are in this file: ESAG.zip.

Once you have unzipped the file, open MATLAB, change directory in MATLAB to the newly created ESAG directory, then the following will work.

Loading and plotting data

First to load some data:

Y = csvread('tasmanianData.csv');

The variable Y is then a 3-by-33 array, the columns being 33 different spherical (unit-vector) observations. They can be plotted with

figure; plotSphere();
plotPointsOnSphere(Y);
view(230,35);

Fitting an ESAG distribution by MLE

An ESAG distribution can be fitted to these data with

parHat = ESAGmle(Y);

which returns parHat, the MLE of the parameter vector, in the form [mu1, mu2, mu3, gamma1, gamma2]. The fitted mean direction is the directional component of mu, i.e.,

meanDir = parHat(1:3)/norm(parHat(1:3));

which we can add to the plot (as a point in yellow) with

plotPointsOnSphere(meanDir,'y');

The ESAG density can be evaluated with the function ESAGdensity, e.g. as follows at the point [0 1 0],

ESAGdensity([0 1 0]', parHat);

Plotting ESAG contours

ESAG distributions can be visualised by plotting contours of constant probability density. The command

plotContoursESAG(parHat,[0.9 0.5 0.1]);

plots contours for the model fitted above, where 0.9, 0.5 and 0.1 are the coverage levels of the contours; e.g., 0.9 is the contour that contains 90% of the probability. Note that this function computes the contour values using Monte Carlo, so it typically takes a few seconds to run.

Simulating from ESAG

One of the big advantages of ESAG is that it is very quick to simulate from. The following simulates n = 100 draws from an ESAG distribution with parameters par

n = 100;
par = [0 0 5 0.8 0];
X = ESAGsim(par,n);

To get a sense of the distribution try visualising these simulated data with plotPointsOnSphere(X); adding contours (e.g. as follows with 0.8 coverage) plotContoursESAG(par,0.8), and repeating with different values of par.

Fitting the isotropic special case, and testing the significance of anisotropy

The "isotropic angular Gaussian" (IAG) distribution, which has circular rather than elliptical contours, is the special case of ESAG with the parameters gam1 = gam2 = 0. It can be fitted easily with ESAGmle (working again with the data loaded above) as follows

parToFit = [1 1 1 0 0];
par0 = [1 1 1 0 0];
[~, fIAG] = ESAGmle(Y,[],par0,parToFit);

in which parToFit specifies which ESAG parameters are to be optimised over - here mu1, mu2, and m3, but not gam1 or gam 2 - and par0 specifies the starting values of the parameters in the optimisation (the important thing is that the gam1 and gam2 are set to zero for IAG; the values of 1 for mu1, mu2, and mu3 here are arbitrary and the optimisation typically converges whatever values are used). The output fIAG is the value of (minus) the maximised log-likelihood for IAG, which we compare with the corresponding value for the general ESAG distribution