--- a/scripts/statistics/base/moment.m	Mon Jan 03 18:36:49 2011 -0800
+++ b/scripts/statistics/base/moment.m	Mon Jan 03 21:23:08 2011 -0800
@@ -18,25 +18,89 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {} moment (@var{x}, @var{p}, @var{opt}, @var{dim})
-## If @var{x} is a vector, compute the @var{p}-th moment of @var{x}.
+## @deftypefn  {Function File} {} moment (@var{x}, @var{p})
+## @deftypefnx {Function File} {} moment (@var{x}, @var{p}, @var{type})
+## @deftypefnx {Function File} {} moment (@var{x}, @var{p}, @var{dim})
+## @deftypefnx {Function File} {} moment (@var{x}, @var{p}, @var{type}, @var{dim})
+## @deftypefnx {Function File} {} moment (@var{x}, @var{p}, @var{dim}, @var{type})
+## Compute the @var{p}-th moment of the vector @var{x} about zero.
+## @tex
+## $$
+## {\rm moment} (x) = { \sum_{i=1}^N {x_i}^p \over N }
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+## moment (x) = 1/N SUM_i x(i)^p
+## @end group
+## @end example
+##
+## @end ifnottex
 ##
 ## If @var{x} is a matrix, return the row vector containing the
 ## @var{p}-th moment of each column.
 ##
-## With the optional string opt, the kind of moment to be computed can
-## be specified.  If opt contains @code{"c"} or @code{"a"}, central
-## and/or absolute moments are returned.  For example,
+## The optional string @var{type} specifies the type of moment to be computed.
+## Valid options are:
+## @table @code
+## @item "c"
+##   Central Moment.  The moment about the mean defined as
+## @tex
+## $$
+## {\sum_{i=1}^N (x_i - \bar{x})^p \over N}
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+## 1/N SUM_i (x(i) - mean(x))^p
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+## @item "a"
+##   Absolute Moment.  The moment about zero ignoring sign defined as
+## @tex
+## $$
+## {\sum_{i=1}^N {\left| x_i \right|}^p \over N}
+## $$
+## @end tex
+## @ifnottex
 ##
 ## @example
-## moment (x, 3, "ac")
+## @group
+## 1/N SUM_i ( abs(x(i)) )^p
+## @end group
 ## @end example
 ##
-## @noindent
-## computes the third central absolute moment of @var{x}.
+## @end ifnottex
+##
+## @item "ac"
+##   Absolute Central Moment.  Defined as
+## @tex
+## $$
+## {\sum_{i=1}^N {\left| x_i - \bar{x} \right|}^p \over N}
+## $$
+## @end tex
+## @ifnottex
 ##
-## If the optional argument @var{dim} is supplied, work along dimension
-## @var{dim}.
+## @example
+## @group
+## 1/N SUM_i ( abs(x(i) - mean(x)) )^p
+## @end group
+## @end example
+##
+## @end ifnottex
+## @end table
+##
+## If the optional argument @var{dim} is given, operate along this dimension.
+##
+## If both @var{type} and @var{dim} are given they may appear in any order.
+## @seealso{var,skewness,kurtosis}
 ## @end deftypefn
 
 ## Can easily be made to work for continuous distributions (using quad)
@@ -51,58 +115,80 @@
     print_usage ();
   endif
 
+  if (!isnumeric(x) || isempty(x) )
+    error ("moment: X must be a non-empty numeric matrix or vector");
+  endif
+
+  if (!(isnumeric(p) && isscalar(p)))
+    error ("moment: P must be a numeric scalar");
+  endif
+
   need_dim = 0;
 
   if (nargin == 2)
-    opt = "";
+    type = "";
     need_dim = 1;
   elseif (nargin == 3)
     if (ischar (opt1))
-      opt = opt1;
+      type = opt1;
       need_dim = 1;
     else
       dim = opt1;
-      opt = "";
+      type = "";
     endif
   elseif (nargin == 4)
     if (ischar (opt1))
-      opt = opt1;
+      type = opt1;
       dim = opt2;
     elseif (ischar (opt2))
-      opt = opt2;
+      type = opt2;
       dim = opt1;
     else
-      error ("moment: expecting opt to be a string");
-    endif
-  else
-    print_usage ();
-  endif
-
-  if (need_dim)
-    t = find (size (x) != 1);
-    if (isempty (t))
-      dim = 1;
-    else
-      dim = t(1);
+      error ("moment: expecting TYPE to be a string");
     endif
   endif
 
+  nd = ndims (x);
   sz = size (x);
-  n = sz (dim);
-
-  if (numel (x) < 1)
-    error ("moment: x must not be empty");
+  if (need_dim)
+    ## Find the first non-singleton dimension.
+    dim = find (sz > 1, 1);
+    if (isempty (dim))
+      dim = 1;
+    endif
+  else
+    if (!(isscalar (dim) && dim == fix (dim)) || 
+        !(1 <= dim && dim <= nd))
+      error ("moment: DIM must be an integer and a valid dimension");
+    endif
   endif
 
-  if any (opt == "c")
+  n = sz(dim);
+
+  if any (type == "c")
     rng = ones (1, length (sz));
     rng(dim) = sz(dim);
     x = x - repmat (sum (x, dim), rng) / n;
   endif
-  if any (opt == "a")
+  if any (type == "a")
     x = abs (x);
   endif
 
   m = sum (x .^ p, dim) / n;
 
 endfunction
+
+
+%% Test input validation
+%!error moment ()
+%!error moment (1)
+%!error moment (1, 2, 3, 4, 5)
+%!error moment ([true true], 2)
+%!error moment (ones(2,0,3), 2)
+%!error moment (1, true)
+%!error moment (1, ones(2,2))
+%!error moment (1, 2, 3, 4)
+%!error moment (1, 2, ones(2,2))
+%!error moment (1, 2, 1.5)
+%!error moment (1, 2, 4)
+