[6463] in s-news-athena
Minimal ellips
daemon@ATHENA.MIT.EDU (Paul Eilers)
Mon Feb 20 04:23:35 1995
From: Paul Eilers <paul@dcmr.nl>
To: s-news@utstat.toronto.edu
Date: Mon, 20 Feb 95 9:45:00 MET
Last week there was a question by Pierre Duchesne about ellipsoidal
hulls. There is a paper by D. M. Titterington that solves this problem
(Applied Statistics 27, 1978, 227-234).
Titterington refers to other papers, by himself and others, but
gives a clear algorithm, for any number of dimensions.
Some time ago I wrote a Matlab script to do the computations. It
follows below. It should not be too difficult to translate it to
S-plus. (In Matlab, X*Y is a matrix product, and X .* Y is an
element by element product. X' is the transpose of X, and
sum(X) gives the sum of the columns of X.)
The results can easily be presented in two dimensions. In three or more
I'm lost.
The speed of the algorithm is very dependent on the data. In some
cases it takes less than 50 iterations to converge, in other cases
over 1000! Titterington mentions the possibly low speed in his paper.
Paul Eilers (paul@dcmr.nl)
DCMR
Schiedam
The Netherlands
===================================================================
% Here is the Matlab script
% It uses the functions "minlips" and "plotlips" defined below
rand('normal');
m = 50;
Y = rand(m, 2);
[T g] = minlips(Y);
plotlips(T, g, 2);
hold;
plot(Y(:, 1), Y(:, 2), 'o');
hold;
end;
function [T , g] = minlips(Y)
% Compute the minimal ellips
% D. M. Titterington, Applied Statistics, 27 (1978) 227 - 234.
[m n] = size(Y);
w = ones(m, 1);
w1 = 0* w;
while max(abs(w - w1)) > 1e-4
w1 = w;
W = (w / sum(w)) * ones(1, n);
R = Y - ones(m, 1) * sum(W .* Y);
S = R' * (W .* R);
D = R .* (R * inv(S));
d = D * ones(n,1);
w = w .* d / n;
disp(max(abs(w - w1)));
end;
T = inv(S);
g = sum(W .* Y);
end;
function plotlips(T, g, c)
% Plot the ellips (u - g)' * T * (u - g) = c
f = linspace(0, 2 * pi, 100)';
U = [sin(f) cos(f)];
D = U .* (U * T);
r = sqrt(c ./ (D * ones(2, 1)));
x = r .* U(:, 1) + g(1);
y = r .* U(:, 2) + g(2);
plot(x, y);
end;
