matthiasm@8: function h=plotgauss2d(mu, Sigma)
matthiasm@8: % PLOTGAUSS2D Plot a 2D Gaussian as an ellipse with optional cross hairs
matthiasm@8: % h=plotgauss2(mu, Sigma)
matthiasm@8: %
matthiasm@8: 
matthiasm@8: h = plotcov2(mu, Sigma);
matthiasm@8: return;
matthiasm@8: 
matthiasm@8: %%%%%%%%%%%%%%%%%%%%%%%%
matthiasm@8: 
matthiasm@8: % PLOTCOV2 - Plots a covariance ellipse with major and minor axes
matthiasm@8: %            for a bivariate Gaussian distribution.
matthiasm@8: %
matthiasm@8: % Usage:
matthiasm@8: %   h = plotcov2(mu, Sigma[, OPTIONS]);
matthiasm@8: % 
matthiasm@8: % Inputs:
matthiasm@8: %   mu    - a 2 x 1 vector giving the mean of the distribution.
matthiasm@8: %   Sigma - a 2 x 2 symmetric positive semi-definite matrix giving
matthiasm@8: %           the covariance of the distribution (or the zero matrix).
matthiasm@8: %
matthiasm@8: % Options:
matthiasm@8: %   'conf'    - a scalar between 0 and 1 giving the confidence
matthiasm@8: %               interval (i.e., the fraction of probability mass to
matthiasm@8: %               be enclosed by the ellipse); default is 0.9.
matthiasm@8: %   'num-pts' - the number of points to be used to plot the
matthiasm@8: %               ellipse; default is 100.
matthiasm@8: %
matthiasm@8: % This function also accepts options for PLOT.
matthiasm@8: %
matthiasm@8: % Outputs:
matthiasm@8: %   h     - a vector of figure handles to the ellipse boundary and
matthiasm@8: %           its major and minor axes
matthiasm@8: %
matthiasm@8: % See also: PLOTCOV3
matthiasm@8: 
matthiasm@8: % Copyright (C) 2002 Mark A. Paskin
matthiasm@8: 
matthiasm@8: function h = plotcov2(mu, Sigma, varargin)
matthiasm@8: 
matthiasm@8: if size(Sigma) ~= [2 2], error('Sigma must be a 2 by 2 matrix'); end
matthiasm@8: if length(mu) ~= 2, error('mu must be a 2 by 1 vector'); end
matthiasm@8: 
matthiasm@8: [p, ...
matthiasm@8:  n, ...
matthiasm@8:  plot_opts] = process_options(varargin, 'conf', 0.9, ...
matthiasm@8: 					'num-pts', 100);
matthiasm@8: h = [];
matthiasm@8: holding = ishold;
matthiasm@8: if (Sigma == zeros(2, 2))
matthiasm@8:   z = mu;
matthiasm@8: else
matthiasm@8:   % Compute the Mahalanobis radius of the ellipsoid that encloses
matthiasm@8:   % the desired probability mass.
matthiasm@8:   k = conf2mahal(p, 2);
matthiasm@8:   % The major and minor axes of the covariance ellipse are given by
matthiasm@8:   % the eigenvectors of the covariance matrix.  Their lengths (for
matthiasm@8:   % the ellipse with unit Mahalanobis radius) are given by the
matthiasm@8:   % square roots of the corresponding eigenvalues.
matthiasm@8:   if (issparse(Sigma))
matthiasm@8:     [V, D] = eigs(Sigma);
matthiasm@8:   else
matthiasm@8:     [V, D] = eig(Sigma);
matthiasm@8:   end
matthiasm@8:   % Compute the points on the surface of the ellipse.
matthiasm@8:   t = linspace(0, 2*pi, n);
matthiasm@8:   u = [cos(t); sin(t)];
matthiasm@8:   w = (k * V * sqrt(D)) * u;
matthiasm@8:   z = repmat(mu, [1 n]) + w;
matthiasm@8:   % Plot the major and minor axes.
matthiasm@8:   L = k * sqrt(diag(D));
matthiasm@8:   h = plot([mu(1); mu(1) + L(1) * V(1, 1)], ...
matthiasm@8: 	   [mu(2); mu(2) + L(1) * V(2, 1)], plot_opts{:});
matthiasm@8:   hold on;
matthiasm@8:   h = [h; plot([mu(1); mu(1) + L(2) * V(1, 2)], ...
matthiasm@8: 	       [mu(2); mu(2) + L(2) * V(2, 2)], plot_opts{:})];
matthiasm@8: end
matthiasm@8: 
matthiasm@8: h = [h; plot(z(1, :), z(2, :), plot_opts{:})];
matthiasm@8: if (~holding) hold off; end
matthiasm@8: 
matthiasm@8: %%%%%%%%%%%%
matthiasm@8: 
matthiasm@8: % CONF2MAHAL - Translates a confidence interval to a Mahalanobis
matthiasm@8: %              distance.  Consider a multivariate Gaussian
matthiasm@8: %              distribution of the form
matthiasm@8: %
matthiasm@8: %   p(x) = 1/sqrt((2 * pi)^d * det(C)) * exp((-1/2) * MD(x, m, inv(C)))
matthiasm@8: %
matthiasm@8: %              where MD(x, m, P) is the Mahalanobis distance from x
matthiasm@8: %              to m under P:
matthiasm@8: %
matthiasm@8: %                 MD(x, m, P) = (x - m) * P * (x - m)'
matthiasm@8: %
matthiasm@8: %              A particular Mahalanobis distance k identifies an
matthiasm@8: %              ellipsoid centered at the mean of the distribution.
matthiasm@8: %              The confidence interval associated with this ellipsoid
matthiasm@8: %              is the probability mass enclosed by it.  Similarly,
matthiasm@8: %              a particular confidence interval uniquely determines
matthiasm@8: %              an ellipsoid with a fixed Mahalanobis distance.
matthiasm@8: %
matthiasm@8: %              If X is an d dimensional Gaussian-distributed vector,
matthiasm@8: %              then the Mahalanobis distance of X is distributed
matthiasm@8: %              according to the Chi-squared distribution with d
matthiasm@8: %              degrees of freedom.  Thus, the Mahalanobis distance is
matthiasm@8: %              determined by evaluating the inverse cumulative
matthiasm@8: %              distribution function of the chi squared distribution
matthiasm@8: %              up to the confidence value.
matthiasm@8: %
matthiasm@8: % Usage:
matthiasm@8: % 
matthiasm@8: %   m = conf2mahal(c, d);
matthiasm@8: %
matthiasm@8: % Inputs:
matthiasm@8: %
matthiasm@8: %   c    - the confidence interval
matthiasm@8: %   d    - the number of dimensions of the Gaussian distribution
matthiasm@8: %
matthiasm@8: % Outputs:
matthiasm@8: %
matthiasm@8: %   m    - the Mahalanobis radius of the ellipsoid enclosing the
matthiasm@8: %          fraction c of the distribution's probability mass
matthiasm@8: %
matthiasm@8: % See also: MAHAL2CONF
matthiasm@8: 
matthiasm@8: % Copyright (C) 2002 Mark A. Paskin
matthiasm@8: 
matthiasm@8: function m = conf2mahal(c, d)
matthiasm@8: 
matthiasm@8: m = chi2inv(c, d); % matlab stats toolbox