lor, 31 01 2009 kl. 12:33 -0500, skrev John W. Eaton:
On 31-Jan-2009, Soren Hauberg wrote:
| Ok. How about
|
| ellipse (a, shift, level, n, ...)
|
| ? I guess that orders the input arguments by how important they are.
I think the relative order is subjective. So I don't really see a
good reason for changing from the original ordering.
I don't entirely agree with the ordering being subjective, but this is
really no big deal, so I've used the original ordering.
Another option is to keep them in the original ordering but allow an
empty matrix to be used as a placeholder that means "use the default
value". Then you could write
ellipse (a, [], [], shift);
is that an acceptable compromise?
The function is implemented using default arguments, so you can do
ellipse (a, :, :, shift)
It's even documented :-)
Soren
------------------------------------------------------------------------
# HG changeset patch
# User Soren Hauberg <address@hidden>
# Date 1233428405 -3600
# Node ID 8b338522286d8208e5e086effb27375ea5dd6853
# Parent 4385bb503467d6cbd834378dd4023b1f5052b858
scripts/plot/ellipse.m: new function
diff -r 4385bb503467 -r 8b338522286d scripts/ChangeLog
--- a/scripts/ChangeLog Thu Jan 29 18:13:06 2009 -0500
+++ b/scripts/ChangeLog Sat Jan 31 20:00:05 2009 +0100
@@ -1,3 +1,7 @@
+2009-01-31 Soren Hauberg <address@hidden>
+
+ * plot/ellipse.m: New function.
+
2009-01-29 John W. Eaton <address@hidden>
* miscellaneous/fileparts.m: Match all possible file separators.
diff -r 4385bb503467 -r 8b338522286d scripts/plot/Makefile.in
--- a/scripts/plot/Makefile.in Thu Jan 29 18:13:06 2009 -0500
+++ b/scripts/plot/Makefile.in Sat Jan 31 20:00:05 2009 +0100
@@ -98,6 +98,7 @@
cylinder.m \
diffuse.m \
gnuplot_drawnow.m \
+ ellipse.m \
ellipsoid.m \
errorbar.m \
ezcontourf.m \
diff -r 4385bb503467 -r 8b338522286d scripts/plot/ellipse.m
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/ellipse.m Sat Jan 31 20:00:05 2009 +0100
@@ -0,0 +1,276 @@
+## Copyright (C) 2001, James B. Rawlings and John W. Eaton
+##
+## 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 3, 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; see the file COPYING. If not, write to
+## the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
+## MA 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} ellipse (@var{a}, @var{level})
+## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{n})
+## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{n},
@var{shift})
+## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{n},
@var{shift}, @dots{})
+## @deftypefnx{Function File} address@hidden =} ellipse (@dots{})
+## @deftypefnx{Function File} address@hidden, @var{y}, @var{major},
@var{minor}, @
+## @var{bbox}] =} ellipse (@dots{})
+##
+## Plot an ellipse with principal axes given the eigenvectors of a 2 by 2
matrix.
+##
+## The eigenvectors of @var{a} multiplied by the corresponding eigenvalues and
+## @var{level} are used as the principal axes of the plotted ellipse. By
default
+## @var{level} is 1. If it is a vector, several ellipses are plotted.
+##
+## The optional argument @var{shift} controls the center of the ellipse. By
default
+## this is a two-vector of zeros.
+##
+## The optional argument @var{n} controls the number of points used to plot the
+## ellipse. By default 100 points are used.
+##
+## The rest of the arguments are passed to @code{plot}, and can be used to
control
+## the visual style of the plot.
+##
+## Besides the usual style arguments that can be passed to plot, @code{ellipse}
+## also accepts the following strings that adds further graphics to the figure.
+##
+## @table @t
+## @item "majoraxis"
+## Also draw the major axis of the ellipse.
+## @item "minoraxis"
+## Also draw the minor axis of the ellipse.
+## @item "boundingbox"
+## Also draw the bounding box of the ellipse.
+## @end table
+##
+## @noindent
+## Input arguments following any of these strings will be used to control the
+## style of the extra graphics elements.
+##
+## If one output argument is requested, a graphics handle of the plot is
returned.
+##
+## If more than one output argument is requested, no plot is generated. Instead
+## the data used to generate the plot is returned. @var{x} and @var{y} are
+## @var{n}-vectors containing the points on the ellipse. @var{major} and
@var{minor}
+## contains the axes of the ellipse, and @var{bbox} contains the bounding box
of
+## the ellipse.
+##
+## The following example generates non-isotropic normally distributed data, and
+## plots the ellipses with Mahalanobis' distance of 1 and 3 to the center.
+##
+## @example
+## ## Generate data
+## X = randn (100, 2);
+## scale = 2;
+## X (:, 1) *= scale;
+## theta = 0.4;
+## R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
+## X = X * R;
+##
+## ## Estimate mean and covariance matrix
+## mu = mean (X);
+## C = cov (X);
+##
+## ## Plot data and ellipses
+## figure
+## plot (X (:, 1), X (:, 2), 'ro');
+## hold on, ellipse (inv (C), [1, 3], :, mu, "k"); hold off
+## axis equal
+## @end example
+##
+## @noindent
+## Note the use of @code{:} to denote the default value of @var{n}. To add a
+## green major axis of width 2, and a blue bounding box to the figure, replace
+## the call to @code{ellipse} with
+##
+## @example
+## ellipse (inv (C), [1, 3], :, mu, "k", "majoraxis", "g", ...
+## "linewidth", 2, "boundingbox", "b");
+## @end example
+## @seealso{plot}
+## @end deftypefn
+
+function varargout = ellipse (semiaxes, level = 1, n = 100, shift = [0, 0],
varargin)
+ ## Check input
+ if (nargin < 1)
+ error ("ellipse: not enough input arguments");
+ endif
+
+ if (!isnumeric (semiaxes) || !isequal (size (semiaxes), [2, 2]))
+ error ("ellipse: first input argument must be a 2x2 matrix");
+ endif
+
+ if (!isvector (level))
+ error ("ellipse: second input argument must be a vector");
+ endif
+
+ if (!isnumeric (shift) || numel (shift) != 2)
+ error ("ellipse: third input argument must be a two-vector");
+ endif
+ shift = shift (:).';
+
+ if (!isscalar (n) || n < 0 || n != round (n))
+ error ("ellipse: fourth input argument must be a positive integer");
+ endif
+
+ ## Possibly override default options
+ available_options = {"minoraxis", "majoraxis", "boundingbox"};
+ num_options = length (available_options);
+ option_idx = zeros (1, num_options);
+ for idx = 1:length (varargin)
+ for k = 1:num_options
+ opt = available_options {k};
+ if (strcmpi (opt, varargin {idx}))
+ option_idx (k) = idx;
+ endif
+ endfor
+ endfor
+
+ [option_idx, idx] = sort (option_idx);
+ available_options = available_options (idx);
+ option_idx = [option_idx, length(varargin)+1];
+
+ ellipse_opt = minoraxis_opt = majoraxis_opt = boundingbox_opt = {};
+
+ last_idx = option_idx (find (option_idx > 0, 1));
+ ellipse_opt = varargin (1:last_idx-1);
+ for k = 1:num_options
+ opt = available_options {k};
+ if (option_idx (k) > 0)
+ val = varargin (last_idx+1:option_idx (k+1)-1);
+ last_idx = option_idx (k+1);
+ eval (sprintf ("%s_opt = val;", opt));
+ endif
+ endfor
+
+ ## We're ready to do the actual work
+ [v, l] = eig (semiaxes);
+
+ dl = diag (l);
+ if (any (imag (dl)) || any (dl <= 0))
+ error ("ellipse: first input argument must be positive definite");
+ endif
+
+ ## Generate contour data.
+ a = 1 / sqrt (dl (1));
+ b = 1 / sqrt (dl (2));
+
+ a = a * level (:);
+ b = b * level (:);
+
+ t = linspace (0, 2*pi, n);
+
+ xt = a * cos (t); xt = xt.';
+ yt = b * sin (t); yt = yt.';
+
+ ## Rotate the contours.
+ ra = atan2 (v (2, 1), v (1, 1));
+
+ cos_ra = cos (ra);
+ sin_ra = sin (ra);
+
+ x = xt * cos_ra - yt * sin_ra + shift (1);
+ y = xt * sin_ra + yt * cos_ra + shift (2);
+
+ ## Endpoints of the major and minor axes.
+ [mx, ii] = max (level);
+
+ minor = (v * diag ([a(ii), b(ii)])).';
+ major = minor;
+
+ major (2, :) = -major (1,:);
+ minor (1, :) = -minor (2,:);
+
+ t = [1; 1] * shift;
+
+ major = major + t;
+ minor = minor + t;
+
+ ## Bounding box for the ellipse using magic formula.
+ ainv = inv (semiaxes);
+ xbox = level (ii) * sqrt (ainv (1,1));
+ ybox = level(ii) * sqrt (ainv (2,2));
+
+ bbox = [xbox ybox; xbox -ybox; -xbox -ybox; -xbox ybox; xbox ybox];
+
+ t = [1; 1; 1; 1; 1] * shift;
+ bbox = bbox + t;
+
+ ## Should we plot the ellipse, or return it?
+ if (nargout <= 1)
+ plot_args = {x, y, ellipse_opt{:}};
+
+ if (!isempty (minoraxis_opt))
+ plot_args = {plot_args{:}, minor(:,1), minor(:,2), minoraxis_opt{:}};
+ endif
+
+ if (!isempty (majoraxis_opt))
+ plot_args = {plot_args{:}, major(:,1), major(:,2), majoraxis_opt{:}};
+ endif
+
+ if (!isempty (boundingbox_opt))
+ plot_args = {plot_args{:}, bbox(:,1), bbox(:,2), boundingbox_opt{:}};
+ endif
+ handle = plot (plot_args {:});
+
+ if (nargout == 1)
+ varargout {1} = handle;
+ endif
+ else
+ varargout {1} = x;
+ varargout {2} = y;
+ varargout {3} = major;
+ varargout {4} = minor;
+ varargout {5} = bbox;
+ endif
+endfunction
+
+%!demo
+%! ## Generate data
+%! X = randn (100, 2);
+%! scale = 2;
+%! X (:, 1) *= scale;
+%! theta = 0.4;
+%! R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
+%! X = X * R;
+%!
+%! ## Estimate mean and covariance matrix
+%! mu = mean (X);
+%! C = cov (X);
+%!
+%! ## Plot data and ellipses
+%! figure
+%! plot (X (:, 1), X (:, 2), 'ro', 'linewidth', 3);
+%! hold on, ellipse (inv (C), [1, 3], :, mu, "k"); hold off
+%! axis equal
+
+%!demo
+%! ## Generate data
+%! X = randn (100, 2);
+%! scale = 2;
+%! X (:, 1) *= scale;
+%! theta = 0.4;
+%! R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
+%! X = X * R;
+%!
+%! ## Estimate mean and covariance matrix
+%! mu = mean (X);
+%! C = cov (X);
+%!
+%! ## Plot data and ellipses. This time also with a major axis and a bounding
box.
+%! figure
+%! plot (X (:, 1), X (:, 2), 'ro', 'linewidth', 3);
+%! hold on
+%! ellipse (inv (C), [1, 3], :, mu, "k", "majoraxis", "g", ...
+%! "linewidth", 2, "boundingbox", "b");
+%! hold off
+%! axis equal
+