wolffd@0: % PLOTCOV2 - Plots a covariance ellipse with major and minor axes
wolffd@0: %            for a bivariate Gaussian distribution.
wolffd@0: %
wolffd@0: % Usage:
wolffd@0: %   h = plotcov2(mu, Sigma[, OPTIONS]);
wolffd@0: % 
wolffd@0: % Inputs:
wolffd@0: %   mu    - a 2 x 1 vector giving the mean of the distribution.
wolffd@0: %   Sigma - a 2 x 2 symmetric positive semi-definite matrix giving
wolffd@0: %           the covariance of the distribution (or the zero matrix).
wolffd@0: %
wolffd@0: % Options:
wolffd@0: %   'conf'    - a scalar between 0 and 1 giving the confidence
wolffd@0: %               interval (i.e., the fraction of probability mass to
wolffd@0: %               be enclosed by the ellipse); default is 0.9.
wolffd@0: %   'num-pts' - the number of points to be used to plot the
wolffd@0: %               ellipse; default is 100.
wolffd@0: %
wolffd@0: % This function also accepts options for PLOT.
wolffd@0: %
wolffd@0: % Outputs:
wolffd@0: %   h     - a vector of figure handles to the ellipse boundary and
wolffd@0: %           its major and minor axes
wolffd@0: %
wolffd@0: % See also: PLOTCOV3
wolffd@0: 
wolffd@0: % Copyright (C) 2002 Mark A. Paskin
wolffd@0: %
wolffd@0: % This program is free software; you can redistribute it and/or modify
wolffd@0: % it under the terms of the GNU General Public License as published by
wolffd@0: % the Free Software Foundation; either version 2 of the License, or
wolffd@0: % (at your option) any later version.
wolffd@0: %
wolffd@0: % This program is distributed in the hope that it will be useful, but
wolffd@0: % WITHOUT ANY WARRANTY; without even the implied warranty of
wolffd@0: % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
wolffd@0: % General Public License for more details.
wolffd@0: %
wolffd@0: % You should have received a copy of the GNU General Public License
wolffd@0: % along with this program; if not, write to the Free Software
wolffd@0: % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
wolffd@0: % USA.
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: 
wolffd@0: function h = plotcov2(mu, Sigma, varargin)
wolffd@0: 
wolffd@0: if size(Sigma) ~= [2 2], error('Sigma must be a 2 by 2 matrix'); end
wolffd@0: if length(mu) ~= 2, error('mu must be a 2 by 1 vector'); end
wolffd@0: 
wolffd@0: [p, ...
wolffd@0:  n, ...
wolffd@0:  plot_opts] = process_options(varargin, 'conf', 0.9, ...
wolffd@0: 					'num-pts', 100);
wolffd@0: h = [];
wolffd@0: holding = ishold;
wolffd@0: if (Sigma == zeros(2, 2))
wolffd@0:   z = mu;
wolffd@0: else
wolffd@0:   % Compute the Mahalanobis radius of the ellipsoid that encloses
wolffd@0:   % the desired probability mass.
wolffd@0:   k = conf2mahal(p, 2);
wolffd@0:   % The major and minor axes of the covariance ellipse are given by
wolffd@0:   % the eigenvectors of the covariance matrix.  Their lengths (for
wolffd@0:   % the ellipse with unit Mahalanobis radius) are given by the
wolffd@0:   % square roots of the corresponding eigenvalues.
wolffd@0:   if (issparse(Sigma))
wolffd@0:     [V, D] = eigs(Sigma);
wolffd@0:   else
wolffd@0:     [V, D] = eig(Sigma);
wolffd@0:   end
wolffd@0:   % Compute the points on the surface of the ellipse.
wolffd@0:   t = linspace(0, 2*pi, n);
wolffd@0:   u = [cos(t); sin(t)];
wolffd@0:   w = (k * V * sqrt(D)) * u;
wolffd@0:   z = repmat(mu, [1 n]) + w;
wolffd@0:   % Plot the major and minor axes.
wolffd@0:   L = k * sqrt(diag(D));
wolffd@0:   h = plot([mu(1); mu(1) + L(1) * V(1, 1)], ...
wolffd@0: 	   [mu(2); mu(2) + L(1) * V(2, 1)], plot_opts{:});
wolffd@0:   hold on;
wolffd@0:   h = [h; plot([mu(1); mu(1) + L(2) * V(1, 2)], ...
wolffd@0: 	       [mu(2); mu(2) + L(2) * V(2, 2)], plot_opts{:})];
wolffd@0: end
wolffd@0: 
wolffd@0: h = [h; plot(z(1, :), z(2, :), plot_opts{:})];
wolffd@0: if (~holding) hold off; end