Mercurial > hg > camir-aes2014
comparison toolboxes/FullBNT-1.0.7/KPMtools/plotcov3.m @ 0:e9a9cd732c1e tip
first hg version after svn
| author | wolffd |
|---|---|
| date | Tue, 10 Feb 2015 15:05:51 +0000 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:e9a9cd732c1e |
|---|---|
| 1 % PLOTCOV3 - Plots a covariance ellipsoid with axes for a trivariate | |
| 2 % Gaussian distribution. | |
| 3 % | |
| 4 % Usage: | |
| 5 % [h, s] = plotcov3(mu, Sigma[, OPTIONS]); | |
| 6 % | |
| 7 % Inputs: | |
| 8 % mu - a 3 x 1 vector giving the mean of the distribution. | |
| 9 % Sigma - a 3 x 3 symmetric positive semi-definite matrix giving | |
| 10 % the covariance of the distribution (or the zero matrix). | |
| 11 % | |
| 12 % Options: | |
| 13 % 'conf' - a scalar between 0 and 1 giving the confidence | |
| 14 % interval (i.e., the fraction of probability mass to | |
| 15 % be enclosed by the ellipse); default is 0.9. | |
| 16 % 'num-pts' - if the value supplied is n, then (n + 1)^2 points | |
| 17 % to be used to plot the ellipse; default is 20. | |
| 18 % 'plot-opts' - a cell vector of arguments to be handed to PLOT3 | |
| 19 % to contol the appearance of the axes, e.g., | |
| 20 % {'Color', 'g', 'LineWidth', 1}; the default is {} | |
| 21 % 'surf-opts' - a cell vector of arguments to be handed to SURF | |
| 22 % to contol the appearance of the ellipsoid | |
| 23 % surface; a nice possibility that yields | |
| 24 % transparency is: {'EdgeAlpha', 0, 'FaceAlpha', | |
| 25 % 0.1, 'FaceColor', 'g'}; the default is {} | |
| 26 % | |
| 27 % Outputs: | |
| 28 % h - a vector of handles on the axis lines | |
| 29 % s - a handle on the ellipsoid surface object | |
| 30 % | |
| 31 % See also: PLOTCOV2 | |
| 32 | |
| 33 % Copyright (C) 2002 Mark A. Paskin | |
| 34 % | |
| 35 % This program is free software; you can redistribute it and/or modify | |
| 36 % it under the terms of the GNU General Public License as published by | |
| 37 % the Free Software Foundation; either version 2 of the License, or | |
| 38 % (at your option) any later version. | |
| 39 % | |
| 40 % This program is distributed in the hope that it will be useful, but | |
| 41 % WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 42 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 43 % General Public License for more details. | |
| 44 % | |
| 45 % You should have received a copy of the GNU General Public License | |
| 46 % along with this program; if not, write to the Free Software | |
| 47 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 | |
| 48 % USA. | |
| 49 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 50 | |
| 51 function [h, s] = plotcov3(mu, Sigma, varargin) | |
| 52 | |
| 53 if size(Sigma) ~= [3 3], error('Sigma must be a 3 by 3 matrix'); end | |
| 54 if length(mu) ~= 3, error('mu must be a 3 by 1 vector'); end | |
| 55 | |
| 56 [p, ... | |
| 57 n, ... | |
| 58 plot_opts, ... | |
| 59 surf_opts] = process_options(varargin, 'conf', 0.9, ... | |
| 60 'num-pts', 20, ... | |
| 61 'plot-opts', {}, ... | |
| 62 'surf-opts', {}); | |
| 63 h = []; | |
| 64 holding = ishold; | |
| 65 if (Sigma == zeros(3, 3)) | |
| 66 z = mu; | |
| 67 else | |
| 68 % Compute the Mahalanobis radius of the ellipsoid that encloses | |
| 69 % the desired probability mass. | |
| 70 k = conf2mahal(p, 3); | |
| 71 % The axes of the covariance ellipse are given by the eigenvectors of | |
| 72 % the covariance matrix. Their lengths (for the ellipse with unit | |
| 73 % Mahalanobis radius) are given by the square roots of the | |
| 74 % corresponding eigenvalues. | |
| 75 if (issparse(Sigma)) | |
| 76 [V, D] = eigs(Sigma); | |
| 77 else | |
| 78 [V, D] = eig(Sigma); | |
| 79 end | |
| 80 if (any(diag(D) < 0)) | |
| 81 error('Invalid covariance matrix: not positive semi-definite.'); | |
| 82 end | |
| 83 % Compute the points on the surface of the ellipsoid. | |
| 84 t = linspace(0, 2*pi, n); | |
| 85 [X, Y, Z] = sphere(n); | |
| 86 u = [X(:)'; Y(:)'; Z(:)']; | |
| 87 w = (k * V * sqrt(D)) * u; | |
| 88 z = repmat(mu(:), [1 (n + 1)^2]) + w; | |
| 89 | |
| 90 % Plot the axes. | |
| 91 L = k * sqrt(diag(D)); | |
| 92 h = plot3([mu(1); mu(1) + L(1) * V(1, 1)], ... | |
| 93 [mu(2); mu(2) + L(1) * V(2, 1)], ... | |
| 94 [mu(3); mu(3) + L(1) * V(3, 1)], plot_opts{:}); | |
| 95 hold on; | |
| 96 h = [h; plot3([mu(1); mu(1) + L(2) * V(1, 2)], ... | |
| 97 [mu(2); mu(2) + L(2) * V(2, 2)], ... | |
| 98 [mu(3); mu(3) + L(2) * V(3, 2)], plot_opts{:})]; | |
| 99 h = [h; plot3([mu(1); mu(1) + L(3) * V(1, 3)], ... | |
| 100 [mu(2); mu(2) + L(3) * V(2, 3)], ... | |
| 101 [mu(3); mu(3) + L(3) * V(3, 3)], plot_opts{:})]; | |
| 102 end | |
| 103 | |
| 104 s = surf(reshape(z(1, :), [(n + 1) (n + 1)]), ... | |
| 105 reshape(z(2, :), [(n + 1) (n + 1)]), ... | |
| 106 reshape(z(3, :), [(n + 1) (n + 1)]), ... | |
| 107 surf_opts{:}); | |
| 108 | |
| 109 if (~holding) hold off; end |