Mercurial > hg > camir-aes2014
diff toolboxes/FullBNT-1.0.7/KPMtools/plotcov2New.m @ 0:e9a9cd732c1e tip
first hg version after svn
author | wolffd |
---|---|
date | Tue, 10 Feb 2015 15:05:51 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/FullBNT-1.0.7/KPMtools/plotcov2New.m Tue Feb 10 15:05:51 2015 +0000 @@ -0,0 +1,119 @@ +% PLOTCOV2 - Plots a covariance ellipsoid with axes for a bivariate +% Gaussian distribution. +% +% Usage: +% [h, s] = plotcov2(mu, Sigma[, OPTIONS]); +% +% Inputs: +% mu - a 2 x 1 vector giving the mean of the distribution. +% Sigma - a 2 x 2 symmetric positive semi-definite matrix giving +% the covariance of the distribution (or the zero matrix). +% +% Options: +% 'conf' - a scalar between 0 and 1 giving the confidence +% interval (i.e., the fraction of probability mass to +% be enclosed by the ellipse); default is 0.9. +% 'num-pts' - if the value supplied is n, then (n + 1)^2 points +% to be used to plot the ellipse; default is 20. +% 'label' - if non-empty, a string that will label the +% ellipsoid (default: []) +% 'plot-axes' - a 0/1 flag indicating if the ellipsoid's axes +% should be plotted (default: 1) +% 'plot-opts' - a cell vector of arguments to be handed to PLOT3 +% to contol the appearance of the axes, e.g., +% {'Color', 'g', 'LineWidth', 1}; the default is {} +% 'fill-color' - a color specifier; is this is not [], the +% covariance ellipse is filled with this color +% (default: []) +% +% Outputs: +% h - a vector of handles on the axis lines +% +% See also: PLOTCOV3 + +% Copyright (C) 2002 Mark A. Paskin +% +% This program is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +% General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +% USA. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [h, s] = plotcov2New(mu, Sigma, varargin) + +h = []; +s = []; + +if size(Sigma) ~= [2 2], error('Sigma must be a 2 by 2 matrix'); end +if length(mu) ~= 2, error('mu must be a 2 by 1 vector'); end + +Sigma = checkpsd(Sigma); + +[p, ... + n, ... + label, ... + plot_axes, ... + plot_opts, ... + fill_color] = process_options(varargin, 'conf', 0.9, ... + 'num-pts', 20, ... + 'label', [], ... + 'plot-axes', 1, ... + 'plot-opts', {}, ... + 'fill-color', []); +holding = ishold; +% Compute the Mahalanobis radius of the ellipsoid that encloses +% the desired probability mass. +k = conf2mahal(p, 2); +% Scale the covariance matrix so the confidence region has unit +% Mahalanobis distance. +Sigma = Sigma * k; +% The axes of the covariance ellipse are given by the eigenvectors of +% the covariance matrix. Their lengths (for the ellipse with unit +% Mahalanobis radius) are given by the square roots of the +% corresponding eigenvalues. +[V, D] = eig(full(Sigma)); +V = real(V); +D = real(D); +D = abs(D); + +% Compute the points on the boundary of the ellipsoid. +t = linspace(0, 2*pi, n); +u = [cos(t(:))'; sin(t(:))']; +w = (V * sqrt(D)) * u; +z = repmat(mu(:), [1 n]) + w; +h = [h; plot(z(1, :), z(2, :), plot_opts{:})]; +if (~isempty(fill_color)) + s = patch(z(1, :), z(2, :), fill_color); +end + +% Plot the axes. +if (plot_axes) + hold on; + L = sqrt(diag(D)); + h = plot([mu(1); mu(1) + L(1) * V(1, 1)], ... + [mu(2); mu(2) + L(1) * V(2, 1)], plot_opts{:}); + h = [h; plot([mu(1); mu(1) + L(2) * V(1, 2)], ... + [mu(2); mu(2) + L(2) * V(2, 2)], plot_opts{:})]; +end + + +if (~isempty(label)) + th = text(mu(1), mu(2), label); + set(th, 'FontSize', 18); + set(th, 'FontName', 'Times'); + set(th, 'FontWeight', 'bold'); + set(th, 'FontAngle', 'italic'); + set(th, 'HorizontalAlignment', 'center'); +end + +if (~holding & plot_axes) hold off; end