# HG changeset patch
# User Julien Bect <julien.bect@supelec.fr>
# Date 1380657077 -7200
#      mar. oct. 01 21:51:17 2013 +0200
# Node ID 0da6ab62aad2d026276024c9f067a2f5d7c0ec87
# Parent  5eae41c8c07f54eb26c4f814bf61de287ff5e5aa
skewness.m: Improve compatibility with Matlab's skewness function

 * skewness.m: Add a 'flag' input argument to select either the sample skewness
   (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 5eae41c8c07f -r 0da6ab62aad2 scripts/statistics/base/skewness.m
--- a/scripts/statistics/base/skewness.m	mar. oct. 01 18:02:23 2013 -0700
+++ b/scripts/statistics/base/skewness.m	mar. oct. 01 21:51:17 2013 +0200
@@ -1,3 +1,4 @@
+## Copyright (C) 2013 Julien Bect
 ## Copyright (C) 1996-2012 John W. Eaton
 ##
 ## This file is part of Octave.
@@ -18,26 +19,34 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {} skewness (@var{x})
-## @deftypefnx {Function File} {} skewness (@var{x}, @var{dim})
-## Compute the skewness of the elements of the vector @var{x}.
+## @deftypefnx {Function File} {} skewness (@var{x}, @var{flag})
+## @deftypefnx {Function File} {} skewness (@var{x}, @var{flag}, @var{dim})
+## Compute the sample skewness of the elements of the vector @var{x}:
 ## @tex
 ## $$
-## {\rm skewness} (x) = {1\over N \sigma^3} \sum_{i=1}^N (x_i-\bar{x})^3
+## {\rm skewness} (x) = {1\over n} \sum_{i=1}^n (x_i-\bar{x})^3 /over \sigma^3
 ## $$
-## where $\bar{x}$ is the mean value of $x$.
+## where $n$ is the length of $x$, $\bar{x}$ its sample mean and $\sigma^2$ its
+## sample variance (biased estimate of the population variance).
 ## @end tex
 ## @ifnottex
 ##
 ## @example
-## skewness (x) = 1/N std(x)^(-3) sum ((x - mean(x)).^3)
+## skewness (@var{X}) = mean ((@var{x} - mean(@var{x})).^3) ./ std(@var{x})^3.
 ## @end example
 ##
 ## @end ifnottex
 ##
 ## @noindent
-## If @var{x} is a matrix, return the skewness along the
-## first non-singleton dimension of the matrix.  If the optional
+## If @var{flag} is equal to 0, return the adjusted skewness coefficient instead
+## of the sample skewness, which is computed using bias-corrected versions of
+## the second and third central moments. If @var{flag} is omitted, empty or
+## equal to 1, return the sample skewness defined above.
+##
+## If @var{x} is a matrix, or more generally a multidimensional array, return
+## the skewness along the first non-singleton dimension.  If the optional
 ## @var{dim} argument is given, operate along this dimension.
+##
 ## @seealso{var, kurtosis, moment}
 ## @end deftypefn
 
@@ -45,57 +54,84 @@
 ## Created: 29 July 1994
 ## Adapted-By: jwe
 
-function retval = skewness (x, dim)
+function y = skewness (x, flag, dim)
 
-  if (nargin != 1 && nargin != 2)
+  if (nargin < 1) || (nargin > 3)
     print_usage ();
   endif
 
-  if (! (isnumeric (x) || islogical (x)))
+  if (! isnumeric (x)) || islogical (x)
     error ("skewness: 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
+    flag = double (flag);
+    if (~ isequal (flag, 0)) && (~ isequal (flag, 1))
+      error ("skewness: FLAG must be 0 or 1");
+    end
+  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 ("skewness: 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 .^ 3, dim);
-  retval(idx) = x(idx) ./ (n * s(idx) .^ 3);
+
+  ## In the following sequence of computations, if x is of class single,
+  ## then the result y is also of class single
+  x = center (x, dim);
+  s = std (x, flag, dim);  # use the biased estimator of the variance
+  y = sum (x .^ 3, dim);
+  y = y ./ (n * s .^ 3);
+  y(s <= 0) = nan;
+
+  ## Apply bias correction to the third central sample moment
+  if flag == 0
+    if n > 2
+      y = y * (n * n ) / ((n - 1) * (n - 2));
+    else
+      y(:) = nan;
+    end
+  endif
 
 endfunction
 
 
-%!assert (skewness ([-1,0,1]), 0)
-%!assert (skewness ([-2,0,1]) < 0)
-%!assert (skewness ([-1,0,2]) > 0)
-%!assert (skewness ([-3,0,1]) == -1*skewness ([-1,0,3]))
+%!assert (skewness ([-1, 0, 1]), 0)
+%!assert (skewness ([-2, 0, 1]) < 0)
+%!assert (skewness ([-1, 0, 2]) > 0)
+%!assert (skewness ([-3, 0, 1]) == -1 * skewness ([-1, 0, 3]))
+%!assert (skewness (ones (3, 5)), nan (1, 5))
+
 %!test
 %! x = [0; 0; 0; 1];
 %! y = [x, 2*x];
-%! assert(all (abs (skewness (y) - [0.75, 0.75]) < sqrt (eps)));
+%! assert (skewness (y),  1.154700538379251 * [1 1], 1e-13)
 
-%!assert (skewness (single (1)), single (0))
+%!assert (skewness ([1:5 10; 1:5 10],  0, 2), 1.439590274527954 * [1; 1], 1e-15)
+%!assert (skewness ([1:5 10; 1:5 10],  1, 2), 1.051328089232020 * [1; 1], 1e-15)
+%!assert (skewness ([1:5 10; 1:5 10], [], 2), 1.051328089232020 * [1; 1], 1e-15)
 
-%% Test input validation
+## Test behaviour on single input
+%! assert (skewness (single ([1:5 10])), single (1.0513283))
+%! assert (skewness (single ([1 2]), 0), single (nan))
+
+## Test input validation
 %!error skewness ()
 %!error skewness (1, 2, 3)
 %!error skewness (['A'; 'B'])
 %!error skewness (1, ones (2,2))
 %!error skewness (1, 1.5)
-%!error skewness (1, 0)
+%!error skewness (1, [], 0)
 %!error skewness (1, 3)
 
