annotate _FullBNT/KPMtools/plotgauss2d.m @ 9:4ea6619cb3f5 tip

removed log files
author matthiasm
date Fri, 11 Apr 2014 15:55:11 +0100
parents b5b38998ef3b
children
rev   line source
matthiasm@8 1 function h=plotgauss2d(mu, Sigma)
matthiasm@8 2 % PLOTGAUSS2D Plot a 2D Gaussian as an ellipse with optional cross hairs
matthiasm@8 3 % h=plotgauss2(mu, Sigma)
matthiasm@8 4 %
matthiasm@8 5
matthiasm@8 6 h = plotcov2(mu, Sigma);
matthiasm@8 7 return;
matthiasm@8 8
matthiasm@8 9 %%%%%%%%%%%%%%%%%%%%%%%%
matthiasm@8 10
matthiasm@8 11 % PLOTCOV2 - Plots a covariance ellipse with major and minor axes
matthiasm@8 12 % for a bivariate Gaussian distribution.
matthiasm@8 13 %
matthiasm@8 14 % Usage:
matthiasm@8 15 % h = plotcov2(mu, Sigma[, OPTIONS]);
matthiasm@8 16 %
matthiasm@8 17 % Inputs:
matthiasm@8 18 % mu - a 2 x 1 vector giving the mean of the distribution.
matthiasm@8 19 % Sigma - a 2 x 2 symmetric positive semi-definite matrix giving
matthiasm@8 20 % the covariance of the distribution (or the zero matrix).
matthiasm@8 21 %
matthiasm@8 22 % Options:
matthiasm@8 23 % 'conf' - a scalar between 0 and 1 giving the confidence
matthiasm@8 24 % interval (i.e., the fraction of probability mass to
matthiasm@8 25 % be enclosed by the ellipse); default is 0.9.
matthiasm@8 26 % 'num-pts' - the number of points to be used to plot the
matthiasm@8 27 % ellipse; default is 100.
matthiasm@8 28 %
matthiasm@8 29 % This function also accepts options for PLOT.
matthiasm@8 30 %
matthiasm@8 31 % Outputs:
matthiasm@8 32 % h - a vector of figure handles to the ellipse boundary and
matthiasm@8 33 % its major and minor axes
matthiasm@8 34 %
matthiasm@8 35 % See also: PLOTCOV3
matthiasm@8 36
matthiasm@8 37 % Copyright (C) 2002 Mark A. Paskin
matthiasm@8 38
matthiasm@8 39 function h = plotcov2(mu, Sigma, varargin)
matthiasm@8 40
matthiasm@8 41 if size(Sigma) ~= [2 2], error('Sigma must be a 2 by 2 matrix'); end
matthiasm@8 42 if length(mu) ~= 2, error('mu must be a 2 by 1 vector'); end
matthiasm@8 43
matthiasm@8 44 [p, ...
matthiasm@8 45 n, ...
matthiasm@8 46 plot_opts] = process_options(varargin, 'conf', 0.9, ...
matthiasm@8 47 'num-pts', 100);
matthiasm@8 48 h = [];
matthiasm@8 49 holding = ishold;
matthiasm@8 50 if (Sigma == zeros(2, 2))
matthiasm@8 51 z = mu;
matthiasm@8 52 else
matthiasm@8 53 % Compute the Mahalanobis radius of the ellipsoid that encloses
matthiasm@8 54 % the desired probability mass.
matthiasm@8 55 k = conf2mahal(p, 2);
matthiasm@8 56 % The major and minor axes of the covariance ellipse are given by
matthiasm@8 57 % the eigenvectors of the covariance matrix. Their lengths (for
matthiasm@8 58 % the ellipse with unit Mahalanobis radius) are given by the
matthiasm@8 59 % square roots of the corresponding eigenvalues.
matthiasm@8 60 if (issparse(Sigma))
matthiasm@8 61 [V, D] = eigs(Sigma);
matthiasm@8 62 else
matthiasm@8 63 [V, D] = eig(Sigma);
matthiasm@8 64 end
matthiasm@8 65 % Compute the points on the surface of the ellipse.
matthiasm@8 66 t = linspace(0, 2*pi, n);
matthiasm@8 67 u = [cos(t); sin(t)];
matthiasm@8 68 w = (k * V * sqrt(D)) * u;
matthiasm@8 69 z = repmat(mu, [1 n]) + w;
matthiasm@8 70 % Plot the major and minor axes.
matthiasm@8 71 L = k * sqrt(diag(D));
matthiasm@8 72 h = plot([mu(1); mu(1) + L(1) * V(1, 1)], ...
matthiasm@8 73 [mu(2); mu(2) + L(1) * V(2, 1)], plot_opts{:});
matthiasm@8 74 hold on;
matthiasm@8 75 h = [h; plot([mu(1); mu(1) + L(2) * V(1, 2)], ...
matthiasm@8 76 [mu(2); mu(2) + L(2) * V(2, 2)], plot_opts{:})];
matthiasm@8 77 end
matthiasm@8 78
matthiasm@8 79 h = [h; plot(z(1, :), z(2, :), plot_opts{:})];
matthiasm@8 80 if (~holding) hold off; end
matthiasm@8 81
matthiasm@8 82 %%%%%%%%%%%%
matthiasm@8 83
matthiasm@8 84 % CONF2MAHAL - Translates a confidence interval to a Mahalanobis
matthiasm@8 85 % distance. Consider a multivariate Gaussian
matthiasm@8 86 % distribution of the form
matthiasm@8 87 %
matthiasm@8 88 % p(x) = 1/sqrt((2 * pi)^d * det(C)) * exp((-1/2) * MD(x, m, inv(C)))
matthiasm@8 89 %
matthiasm@8 90 % where MD(x, m, P) is the Mahalanobis distance from x
matthiasm@8 91 % to m under P:
matthiasm@8 92 %
matthiasm@8 93 % MD(x, m, P) = (x - m) * P * (x - m)'
matthiasm@8 94 %
matthiasm@8 95 % A particular Mahalanobis distance k identifies an
matthiasm@8 96 % ellipsoid centered at the mean of the distribution.
matthiasm@8 97 % The confidence interval associated with this ellipsoid
matthiasm@8 98 % is the probability mass enclosed by it. Similarly,
matthiasm@8 99 % a particular confidence interval uniquely determines
matthiasm@8 100 % an ellipsoid with a fixed Mahalanobis distance.
matthiasm@8 101 %
matthiasm@8 102 % If X is an d dimensional Gaussian-distributed vector,
matthiasm@8 103 % then the Mahalanobis distance of X is distributed
matthiasm@8 104 % according to the Chi-squared distribution with d
matthiasm@8 105 % degrees of freedom. Thus, the Mahalanobis distance is
matthiasm@8 106 % determined by evaluating the inverse cumulative
matthiasm@8 107 % distribution function of the chi squared distribution
matthiasm@8 108 % up to the confidence value.
matthiasm@8 109 %
matthiasm@8 110 % Usage:
matthiasm@8 111 %
matthiasm@8 112 % m = conf2mahal(c, d);
matthiasm@8 113 %
matthiasm@8 114 % Inputs:
matthiasm@8 115 %
matthiasm@8 116 % c - the confidence interval
matthiasm@8 117 % d - the number of dimensions of the Gaussian distribution
matthiasm@8 118 %
matthiasm@8 119 % Outputs:
matthiasm@8 120 %
matthiasm@8 121 % m - the Mahalanobis radius of the ellipsoid enclosing the
matthiasm@8 122 % fraction c of the distribution's probability mass
matthiasm@8 123 %
matthiasm@8 124 % See also: MAHAL2CONF
matthiasm@8 125
matthiasm@8 126 % Copyright (C) 2002 Mark A. Paskin
matthiasm@8 127
matthiasm@8 128 function m = conf2mahal(c, d)
matthiasm@8 129
matthiasm@8 130 m = chi2inv(c, d); % matlab stats toolbox