octave_packages/m/statistics/base/moment.m

   1 ## Copyright (C) 1995-2012 Kurt Hornik
   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} {} moment (@var{x}, @var{p})
  21 ## @deftypefnx {Function File} {} moment (@var{x}, @var{p}, @var{type})
  22 ## @deftypefnx {Function File} {} moment (@var{x}, @var{p}, @var{dim})
  23 ## @deftypefnx {Function File} {} moment (@var{x}, @var{p}, @var{type}, @var{dim})
  24 ## @deftypefnx {Function File} {} moment (@var{x}, @var{p}, @var{dim}, @var{type})
  25 ## Compute the @var{p}-th moment of the vector @var{x} about zero.
  26 ## @tex
  27 ## $$
  28 ## {\rm moment} (x) = { \sum_{i=1}^N {x_i}^p \over N }
  29 ## $$
  30 ## @end tex
  31 ## @ifnottex
  32 ##
  33 ## @example
  34 ## @group
  35 ## moment (x) = 1/N SUM_i x(i)^p
  36 ## @end group
  37 ## @end example
  38 ##
  39 ## @end ifnottex
  40 ##
  41 ## If @var{x} is a matrix, return the row vector containing the
  42 ## @var{p}-th moment of each column.
  43 ##
  44 ## The optional string @var{type} specifies the type of moment to be computed.
  45 ## Valid options are:
  46 ##
  47 ## @table @asis
  48 ## @item "c"
  49 ##   Central Moment.  The moment about the mean defined as
  50 ## @tex
  51 ## $$
  52 ## {\sum_{i=1}^N (x_i - \bar{x})^p \over N}
  53 ## $$
  54 ## @end tex
  55 ## @ifnottex
  56 ##
  57 ## @example
  58 ## @group
  59 ## 1/N SUM_i (x(i) - mean(x))^p
  60 ## @end group
  61 ## @end example
  62 ##
  63 ## @end ifnottex
  64 ##
  65 ## @item "a"
  66 ##   Absolute Moment.  The moment about zero ignoring sign defined as
  67 ## @tex
  68 ## $$
  69 ## {\sum_{i=1}^N {\left| x_i \right|}^p \over N}
  70 ## $$
  71 ## @end tex
  72 ## @ifnottex
  73 ##
  74 ## @example
  75 ## @group
  76 ## 1/N SUM_i ( abs (x(i)) )^p
  77 ## @end group
  78 ## @end example
  79 ##
  80 ## @end ifnottex
  81 ##
  82 ## @item "ac"
  83 ##   Absolute Central Moment.  Defined as
  84 ## @tex
  85 ## $$
  86 ## {\sum_{i=1}^N {\left| x_i - \bar{x} \right|}^p \over N}
  87 ## $$
  88 ## @end tex
  89 ## @ifnottex
  90 ##
  91 ## @example
  92 ## @group
  93 ## 1/N SUM_i ( abs (x(i) - mean(x)) )^p
  94 ## @end group
  95 ## @end example
  96 ##
  97 ## @end ifnottex
  98 ## @end table
  99 ##
 100 ## If the optional argument @var{dim} is given, operate along this dimension.
 101 ##
 102 ## If both @var{type} and @var{dim} are given they may appear in any order.
 103 ## @seealso{var, skewness, kurtosis}
 104 ## @end deftypefn
 105
 106 ## Can easily be made to work for continuous distributions (using quad)
 107 ## as well, but how does the general case work?
 108
 109 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
 110 ## Description: Compute moments
 111
 112 function m = moment (x, p, opt1, opt2)
 113
 114   if (nargin < 2 || nargin > 4)
 115     print_usage ();
 116   endif
 117
 118   if (!(isnumeric (x) || islogical (x)) || isempty (x))
 119     error ("moment: X must be a non-empty numeric matrix or vector");
 120   endif
 121
 122   if (! (isnumeric (p) && isscalar (p)))
 123     error ("moment: P must be a numeric scalar");
 124   endif
 125
 126   need_dim = false;
 127
 128   if (nargin == 2)
 129     type = "";
 130     need_dim = true;
 131   elseif (nargin == 3)
 132     if (ischar (opt1))
 133       type = opt1;
 134       need_dim = true;
 135     else
 136       dim = opt1;
 137       type = "";
 138     endif
 139   elseif (nargin == 4)
 140     if (ischar (opt1))
 141       type = opt1;
 142       dim = opt2;
 143     elseif (ischar (opt2))
 144       type = opt2;
 145       dim = opt1;
 146     else
 147       error ("moment: expecting TYPE to be a string");
 148     endif
 149   endif
 150
 151   nd = ndims (x);
 152   sz = size (x);
 153   if (need_dim)
 154     ## Find the first non-singleton dimension.
 155     (dim = find (sz > 1, 1)) || (dim = 1);
 156   else
 157     if (!(isscalar (dim) && dim == fix (dim)) ||
 158         !(1 <= dim && dim <= nd))
 159       error ("moment: DIM must be an integer and a valid dimension");
 160     endif
 161   endif
 162
 163   n = sz(dim);
 164
 165   if (any (type == "c"))
 166     x = center (x, dim);
 167   endif
 168   if any (type == "a")
 169     x = abs (x);
 170   endif
 171
 172   m = sum (x .^ p, dim) / n;
 173
 174 endfunction
 175
 176
 177 %!test
 178 %! x = rand (10);
 179 %! assert (moment (x,1), mean (x), 1e1*eps);
 180 %! assert (moment (x,2), meansq (x), 1e1*eps);
 181 %! assert (moment (x,1,2), mean (x,2), 1e1*eps);
 182 %! assert (moment (x,1,'c'), mean (center (x)), 1e1*eps);
 183 %! assert (moment (x,1,'a'), mean (abs (x)), 1e1*eps);
 184
 185 %!assert (moment (single([1 2 3]),1), single(2));
 186
 187 %% Test input validation
 188 %!error moment ()
 189 %!error moment (1)
 190 %!error moment (1, 2, 3, 4, 5)
 191 %!error moment (['A'; 'B'], 2)
 192 %!error moment (ones(2,0,3), 2)
 193 %!error moment (1, true)
 194 %!error moment (1, ones(2,2))
 195 %!error moment (1, 2, 3, 4)
 196 %!error moment (1, 2, ones(2,2))
 197 %!error moment (1, 2, 1.5)
 198 %!error moment (1, 2, 4)
 199