octave_packages/m/statistics/base/kurtosis.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} {} kurtosis (@var{x})
  21 ## @deftypefnx {Function File} {} kurtosis (@var{x}, @var{dim})
  22 ## Compute the kurtosis of the elements of the vector @var{x}.
  23 ## @tex
  24 ## $$
  25 ##  {\rm kurtosis} (x) = {1\over N \sigma^4} \sum_{i=1}^N (x_i-\bar{x})^4 - 3
  26 ## $$
  27 ## where $\bar{x}$ is the mean value of $x$.
  28 ## @end tex
  29 ## @ifnottex
  30 ##
  31 ## @example
  32 ## kurtosis (x) = 1/N std(x)^(-4) sum ((x - mean(x)).^4) - 3
  33 ## @end example
  34 ##
  35 ## @end ifnottex
  36 ## If @var{x} is a matrix, return the kurtosis over the
  37 ## first non-singleton dimension of the matrix.  If the optional
  38 ## @var{dim} argument is given, operate along this dimension.
  39 ##
  40 ## Note: The definition of kurtosis above yields a kurtosis of zero for the
  41 ## stdnormal distribution and is sometimes referred to as "excess kurtosis".
  42 ## To calculate kurtosis without the normalization factor of @math{-3} use
  43 ## @code{moment (@var{x}, 4, 'c') / std (@var{x})^4}.
  44 ## @seealso{var, skewness, moment}
  45 ## @end deftypefn
  46
  47 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
  48 ## Created: 29 July 1994
  49 ## Adapted-By: jwe
  50
  51 function retval = kurtosis (x, dim)
  52
  53   if (nargin != 1 && nargin != 2)
  54     print_usage ();
  55   endif
  56
  57   if (! (isnumeric (x) || islogical (x)))
  58     error ("kurtosis: X must be a numeric vector or matrix");
  59   endif
  60
  61   nd = ndims (x);
  62   sz = size (x);
  63   if (nargin != 2)
  64     ## Find the first non-singleton dimension.
  65     (dim = find (sz > 1, 1)) || (dim = 1);
  66   else
  67     if (!(isscalar (dim) && dim == fix (dim))
  68         || !(1 <= dim && dim <= nd))
  69       error ("kurtosis: DIM must be an integer and a valid dimension");
  70     endif
  71   endif
  72
  73   n = sz(dim);
  74   sz(dim) = 1;
  75   x = center (x, dim);  # center also promotes integer to double for next line
  76   retval = zeros (sz, class (x));
  77   s = std (x, [], dim);
  78   idx = find (s > 0);
  79   x = sum (x.^4, dim);
  80   retval(idx) = x(idx) ./ (n * s(idx) .^ 4) - 3;
  81
  82 endfunction
  83
  84
  85 %!test
  86 %! x = [-1; 0; 0; 0; 1];
  87 %! y = [x, 2*x];
  88 %! assert (kurtosis (y), [-1.4, -1.4], sqrt (eps));
  89
  90 %!assert (kurtosis (single(1)), single(0));
  91
  92 %% Test input validation
  93 %!error kurtosis ()
  94 %!error kurtosis (1, 2, 3)
  95 %!error kurtosis (['A'; 'B'])
  96 %!error kurtosis (1, ones(2,2))
  97 %!error kurtosis (1, 1.5)
  98 %!error kurtosis (1, 0)
  99 %!error kurtosis (1, 3)
 100