# HG changeset patch
# User Julien Bect <julien.bect@supelec.fr>
# Date 1382095436 -7200
#      ven. oct. 18 13:23:56 2013 +0200
# Node ID 427c5545ef3fbc99dcb3554b00bf60a1104ac325
# Parent  93e272018df26335dc8fac843fa7212c29431d82
kurtosis.m: Improve compatibility with Matlab's kurtosis function

* kurtosis.m: Add a 'flag' input argument to select either the sample
kurtosis (default) or its "bias corrected" version. Return NaN (instead
of 0) when the standard deviation is zero. Fix documentation. Fix old
tests and add some new ones.

diff -r 93e272018df2 -r 427c5545ef3f scripts/statistics/base/kurtosis.m
--- a/scripts/statistics/base/kurtosis.m	ven. oct. 18 01:36:31 2013 +0200
+++ b/scripts/statistics/base/kurtosis.m	ven. oct. 18 13:23:56 2013 +0200
@@ -1,3 +1,4 @@
+## Copyright (C) 2013 Julien Bect
 ## Copyright (C) 1996-2012 John W. Eaton
 ##
 ## This file is part of Octave.
@@ -52,9 +53,9 @@
 ## Created: 29 July 1994
 ## Adapted-By: jwe
 
-function retval = kurtosis (x, dim)
+function y = kurtosis (x, flag, dim)
 
-  if (nargin != 1 && nargin != 2)
+  if (nargin < 1) || (nargin > 3)
     print_usage ();
   endif
 
@@ -62,26 +63,44 @@
     error ("kurtosis: X must be a numeric vector or matrix");
   endif
 
+  if (nargin < 2 || isempty (flag))
+    flag = 1;  # default: do not use the "bias corrected" version
+  else
+    if ((! isscalar (flag)) || (flag != 0 && flag != 1))
+      error ("kurtosis: FLAG must be 0 or 1");
+    endif
+  endif
+
   nd = ndims (x);
   sz = size (x);
-  if (nargin != 2)
+  if (nargin < 3)
     ## Find the first non-singleton dimension.
     (dim = find (sz > 1, 1)) || (dim = 1);
   else
-    if (!(isscalar (dim) && dim == fix (dim))
-        || !(1 <= dim && dim <= nd))
+    if (! (isscalar (dim) && dim == fix (dim)) || ! (1 <= dim && dim <= nd))
       error ("kurtosis: DIM must be an integer and a valid dimension");
     endif
   endif
 
   n = sz(dim);
   sz(dim) = 1;
-  x = center (x, dim);  # center also promotes integer to double for next line
-  retval = zeros (sz, class (x));
-  s = std (x, [], dim);
-  idx = find (s > 0);
-  x = sum (x.^4, dim);
-  retval(idx) = x(idx) ./ (n * s(idx) .^ 4) - 3;
+
+  x = center (x, dim);   # center also promotes integer, logical to double
+  v = var (x, 1, dim);   # normalize with 1/N
+  y = sum (x .^ 4, dim);
+  idx = (v != 0);
+  y(idx) = y(idx) ./ (n * v(idx) .^ 2);
+  y(! idx) = NaN;
+
+  ## Apply bias correction to the second and fourth central sample moment  
+  if (flag == 0)
+    if (n > 3)
+      C = (n - 1) / ((n - 2) * (n - 3));
+      y(idx) = 3 + C * ((n + 1) * y(idx) - 3 * (n - 1));
+    else
+      y(:) = NaN;
+    endif
+  endif
 
 endfunction
 
@@ -89,16 +108,37 @@
 %!test
 %! x = [-1; 0; 0; 0; 1];
 %! y = [x, 2*x];
-%! assert (kurtosis (y), [-1.4, -1.4], sqrt (eps));
+%! assert (kurtosis (y), [2.5, 2.5], sqrt (eps));
 
-%!assert (kurtosis (single (1)), single (0))
+%!assert (kurtosis ([-3, 0, 1]) == kurtosis ([-1, 0, 3]))
+%!assert (kurtosis (ones (3, 5)), NaN (1, 5))
+
+%!assert (kurtosis ([1:5 10; 1:5 10],  0, 2), 5.437731792528892 * [1; 1], eps)
+%!assert (kurtosis ([1:5 10; 1:5 10],  1, 2), 2.978650900295620 * [1; 1], eps)
+%!assert (kurtosis ([1:5 10; 1:5 10], [], 2), 2.978650900295620 * [1; 1], eps)
+
+## Test behaviour on single input
+%!assert (kurtosis (single ([1:5 10])), single (2.9786513), eps ("single"))
+%!assert (kurtosis (single ([1 2]), 0), single (NaN))
+
+## Verify no "division-by-zero" warnings
+%!test
+%! wstate = warning ("query", "Octave:divide-by-zero");
+%! unwind_protect
+%!   lastwarn ("");  # clear last warning
+%!   kurtosis (1);
+%!   assert (lastwarn (), "");
+%! unwind_protect_cleanup
+%!   warning (wstate, "Octave:divide-by-zero");
+%! end_unwind_protect
 
 %% Test input validation
 %!error kurtosis ()
 %!error kurtosis (1, 2, 3)
-%!error kurtosis (['A'; 'B'])
-%!error kurtosis (1, ones (2,2))
-%!error kurtosis (1, 1.5)
-%!error kurtosis (1, 0)
-%!error kurtosis (1, 3)
-
+%!error <X must be a numeric vector or matrix> kurtosis (['A'; 'B'])
+%!error <FLAG must be 0 or 1> kurtosis (1, 2)
+%!error <FLAG must be 0 or 1> kurtosis (1, [1 0])
+%!error <DIM must be an integer> kurtosis (1, [], ones (2,2))
+%!error <DIM must be an integer> kurtosis (1, [], 1.5)
+%!error <DIM must be .* a valid dimension> kurtosis (1, [], 0)
+%!error <DIM must be .* a valid dimension> kurtosis (1, [], 3)
