--- a/scripts/statistics/base/ols.m
+++ b/scripts/statistics/base/ols.m
@@ -48,9 +48,17 @@
 ##
 ## @table @var
 ## @item beta
-## The OLS estimator for @var{b}, @code{@var{beta} = pinv (@var{x}) *
-## @var{y}}, where @code{pinv (@var{x})} denotes the pseudoinverse of
-## @var{x}.
+## The OLS estimator for @math{b}.
+## @tex
+## $beta$ is calculated directly via $(x^Tx)^{-1} x^T y$ if the matrix $x^Tx$ is
+## of full rank.
+## @end tex
+## @ifnottex
+## @var{beta} is calculated directly via @code{inv (x'*x) * x' * y} if the
+## matrix @code{x'*x} is of full rank.
+## @end ifnottex
+## Otherwise, @code{@var{beta} = pinv (@var{x}) * @var{y}} where 
+## @code{pinv (@var{x})} denotes the pseudoinverse of @var{x}.
 ##
 ## @item sigma
 ## The OLS estimator for the matrix @var{s},
@@ -66,34 +74,63 @@
 ## @item r
 ## The matrix of OLS residuals, @code{@var{r} = @var{y} - @var{x}*@var{beta}}.
 ## @end table
+## @seealso{gls, pinv}
 ## @end deftypefn
 
 ## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at>
 ## Created: May 1993
 ## Adapted-By: jwe
 
-function [BETA, SIGMA, R] = ols (Y, X)
+function [beta, sigma, r] = ols (y, x)
 
   if (nargin != 2)
     print_usage ();
   endif
 
-  [nr, nc] = size (X);
-  [ry, cy] = size (Y);
-  if (nr != ry)
-    error ("ols: incorrect matrix dimensions");
+  if (! (isnumeric (x) && isnumeric (y)))
+    error ("ols: X and Y must be numeric matrices or vectors");
+  endif
+
+  if (ndims (x) != 2 || ndims (y) != 2)
+    error ("ols: X and Y must be 2-D matrices or vectors");
   endif
 
-  Z = X' * X;
-  r = rank (Z);
+  [nr, nc] = size (x);
+  [ry, cy] = size (y);
+  if (nr != ry)
+    error ("ols: number of rows of X and Y must be equal");
+  endif
+
+  z = x' * x;
+  r = rank (z);
 
   if (r == nc)
-    BETA = inv (Z) * X' * Y;
+    beta = inv (z) * x' * y;
   else
-    BETA = pinv (X) * Y;
+    beta = pinv (x) * y;
+  endif
+
+  if (isargout (2) || isargout (3))
+    r = y - x * beta;
+  endif
+  if (isargout (2))
+    sigma = r' * r / (nr - r);
   endif
 
-  R = Y - X * BETA;
-  SIGMA = R' * R / (nr - r);
+endfunction
+
+%!test
+%! x = [1:5]';
+%! y = 3*x + 2;
+%! x = [x, ones(5,1)];
+%! assert (ols(y,x), [3; 2], 50*eps)
 
-endfunction
+%% Test input validation
+%!error ols ();
+%!error ols (1);
+%!error ols (1, 2, 3);
+%!error ols ([true, true], [1, 2]);
+%!error ols ([1, 2], [true, true]);
+%!error ols (ones (2,2,2), ones (2,2));
+%!error ols (ones (2,2), ones (2,2,2));
+%!error ols (ones(1,2), ones(2,2));