octave_packages/m/statistics/base/mahalanobis.m

   1 ## Copyright (C) 1996-2012 John W. Eaton
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} mahalanobis (@var{x}, @var{y})
  21 ## Return the Mahalanobis' D-square distance between the multivariate
  22 ## samples @var{x} and @var{y}, which must have the same number of
  23 ## components (columns), but may have a different number of observations
  24 ## (rows).
  25 ## @end deftypefn
  26
  27 ## Author: Friedrich Leisch <leisch@ci.tuwien.ac.at>
  28 ## Created: July 1993
  29 ## Adapted-By: jwe
  30
  31 function retval = mahalanobis (x, y)
  32
  33   if (nargin != 2)
  34     print_usage ();
  35   endif
  36
  37   if (   ! (isnumeric (x) || islogical (x))
  38       || ! (isnumeric (y) || islogical (y)))
  39     error ("mahalanobis: X and Y must be numeric matrices or vectors");
  40   endif
  41
  42   if (ndims (x) != 2 || ndims (y) != 2)
  43     error ("mahalanobis: X and Y must be 2-D matrices or vectors");
  44   endif
  45
  46   [xr, xc] = size (x);
  47   [yr, yc] = size (y);
  48
  49   if (xc != yc)
  50     error ("mahalanobis: X and Y must have the same number of columns");
  51   endif
  52
  53   if (isinteger (x))
  54     x = double (x);
  55   endif
  56
  57   xm = mean (x);
  58   ym = mean (y);
  59
  60   ## Center data by subtracting means
  61   x = bsxfun (@minus, x, xm);
  62   y = bsxfun (@minus, y, ym);
  63
  64   w = (x' * x + y' * y) / (xr + yr - 2);
  65
  66   winv = inv (w);
  67
  68   retval = (xm - ym) * winv * (xm - ym)';
  69
  70 endfunction
  71
  72
  73 %% Test input validation
  74 %!error mahalanobis ();
  75 %!error mahalanobis (1, 2, 3);
  76 %!error mahalanobis ('A', 'B');
  77 %!error mahalanobis ([1, 2], ['A', 'B']);
  78 %!error mahalanobis (ones (2,2,2));
  79 %!error mahalanobis (ones (2,2), ones (2,2,2));
  80 %!error mahalanobis (ones (2,2), ones (2,3));