--- a/scripts/statistics/base/corr.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/corr.m	Wed Dec 27 23:38:25 2017 -0600
@@ -27,13 +27,13 @@
 ## @var{i}-th variable in @var{x} and the @var{j}-th variable in @var{y}.
 ## @tex
 ## $$
-## {\rm corr}(x,y) = {{\rm cov}(x,y) \over {\rm std}(x) {\rm std}(y)}
+## {\rm corr}(x,y) = {{\rm cov}(x,y) \over {\rm std}(x) \, {\rm std}(y)}
 ## $$
 ## @end tex
 ## @ifnottex
 ##
 ## @example
-## corr (x,y) = cov (x,y) / (std (x) * std (y))
+## corr (@var{x},@var{y}) = cov (@var{x},@var{y}) / (std (@var{x}) * std (@var{y}))
 ## @end example
 ##
 ## @end ifnottex

--- a/scripts/statistics/base/cov.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/cov.m	Wed Dec 27 23:38:25 2017 -0600
@@ -31,14 +31,16 @@
 ## $$
 ## \sigma_{ij} = {1 \over N-1} \sum_{i=1}^N (x_i - \bar{x})(y_i - \bar{y})
 ## $$
-## where $\bar{x}$ and $\bar{y}$ are the mean values of $x$ and $y$.
+## where $\bar{x}$ and $\bar{y}$ are the mean values of @var{x} and @var{y}.
 ## @end tex
 ## @ifnottex
 ##
 ## @example
-## cov (x) = 1/N-1 * SUM_i (x(i) - mean(x)) * (y(i) - mean(y))
+## cov (@var{x}) = 1/(N-1) * SUM_i (@var{x}(i) - mean(@var{x})) * (@var{y}(i) - mean(@var{y}))
 ## @end example
 ##
+## where @math{N} is the length of the @var{x} and @var{y} vectors.
+##
 ## @end ifnottex
 ##
 ## If called with one argument, compute @code{cov (@var{x}, @var{x})}, the

--- a/scripts/statistics/base/gls.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/gls.m	Wed Dec 27 23:38:25 2017 -0600
@@ -18,44 +18,64 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {} {[@var{beta}, @var{v}, @var{r}] =} gls (@var{y}, @var{x}, @var{o})
-## Generalized least squares model.
+## Generalized least squares (GLS) model.
 ##
 ## Perform a generalized least squares estimation for the multivariate model
 ## @tex
-## $y = x b + e$
-## with $\bar{e} = 0$ and cov(vec($e$)) = $(s^2)o$,
+## $@var{y} = @var{x}\,@var{b} + @var{e}$
 ## @end tex
 ## @ifnottex
-## @w{@math{y = x*b + e}} with @math{mean (e) = 0} and
-## @math{cov (vec (e)) = (s^2) o},
+## @w{@math{@var{y} = @var{x}*@var{B} + @var{E}}}
 ## @end ifnottex
 ## where
 ## @tex
-## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix, $b$ is a $k
-## \times p$ matrix, $e$ is a $t \times p$ matrix, and $o$ is a $tp \times
-## tp$ matrix.
+## $@var{y}$ is a $t \times p$ matrix, $@var{x}$ is a $t \times k$ matrix,
+## $@var{b}$ is a $k \times p$ matrix and $@var{e}$ is a $t \times p$ matrix.
 ## @end tex
 ## @ifnottex
-## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by
-## @math{k} matrix, @math{b} is a @math{k} by @math{p} matrix, @math{e}
-## is a @math{t} by @math{p} matrix, and @math{o} is a @math{t*p} by
-## @math{t*p} matrix.
+## @var{y} is a @math{t}-by-@math{p} matrix, @var{x} is a
+## @math{t}-by-@math{k} matrix, @var{b} is a @math{k}-by-@math{p} matrix
+## and @var{e} is a @math{t}-by-@math{p} matrix.
 ## @end ifnottex
 ##
 ## @noindent
-## Each row of @var{y} and @var{x} is an observation and each column a
-## variable.  The return values @var{beta}, @var{v}, and @var{r} are
+## Each row of @var{y} is a @math{p}-variate observation in which each column
+## represents a variable.  Likewise, the rows of @var{x} represent
+## @math{k}-variate observations or possibly designed values.  Furthermore,
+## the collection of observations @var{x} must be of adequate rank, @math{k},
+## otherwise @var{b} cannot be uniquely estimated.
+##
+## The observation errors, @var{e}, are assumed to originate from an
+## underlying @math{p}-variate distribution with zero mean but possibly
+## heteroscedastic observations.  That is, in general,
+## @tex
+## $\bar{@var{e}} = 0$ and cov(vec(@var{e})) = $s^2@var{o}$
+## @end tex
+## @ifnottex
+## @code{@math{mean (@var{e}) = 0}} and
+## @code{@math{cov (vec (@var{e})) = (@math{s}^2)*@var{o}}}
+## @end ifnottex
+## in which @math{s} is a scalar and @var{o} is a
+## @tex
+## @math{t \, p \times t \, p}
+## @end tex
+## @ifnottex
+## @math{t*p}-by-@math{t*p}
+## @end ifnottex
+## matrix.
+##
+## The return values @var{beta}, @var{v}, and @var{r} are
 ## defined as follows.
 ##
 ## @table @var
 ## @item beta
-## The GLS estimator for @math{b}.
+## The GLS estimator for matrix @var{b}.
 ##
 ## @item v
-## The GLS estimator for @math{s^2}.
+## The GLS estimator for scalar @math{s^2}.
 ##
 ## @item r
-## The matrix of GLS residuals, @math{r = y - x*beta}.
+## The matrix of GLS residuals, @math{@var{r} = @var{y} - @var{x}*@var{beta}}.
 ## @end table
 ## @seealso{ols}
 ## @end deftypefn

--- a/scripts/statistics/base/histc.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/histc.m	Wed Dec 27 23:38:25 2017 -0600
@@ -26,8 +26,20 @@
 ## When @var{x} is a vector, the function counts the number of elements of
 ## @var{x} that fall in the histogram bins defined by @var{edges}.  This
 ## must be a vector of monotonically increasing values that define the edges
-## of the histogram bins.  @code{@var{n}(k)} contains the number of elements
-## in @var{x} for which @code{@var{edges}(k) <= @var{x} < @var{edges}(k+1)}.
+## of the histogram bins.
+## @tex
+## $n(k)$
+## @end tex
+## @ifnottex
+## @code{@var{n}(k)}
+## @end ifnottex
+## contains the number of elements in @var{x} for which
+## @tex
+## $@var{edges}(k) <= @var{x} < @var{edges}(k+1)$.
+## @end tex
+## @ifnottex
+## @code{@var{edges}(k) <= @var{x} < @var{edges}(k+1)}.
+## @end ifnottex
 ## The final element of @var{n} contains the number of elements of @var{x}
 ## exactly equal to the last element of @var{edges}.
 ##

--- a/scripts/statistics/base/kendall.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/kendall.m	Wed Dec 27 23:38:25 2017 -0600
@@ -20,23 +20,35 @@
 ## @deftypefn  {} {} kendall (@var{x})
 ## @deftypefnx {} {} kendall (@var{x}, @var{y})
 ## @cindex Kendall's Tau
-## Compute Kendall's @var{tau}.
+## Compute Kendall's
+## @tex
+## $\tau$.
+## @end tex
+## @ifnottex
+## @var{tau}.
+## @end ifnottex
 ##
-## For two data vectors @var{x}, @var{y} of common length @var{n}, Kendall's
-## @var{tau} is the correlation of the signs of all rank differences of
+## For two data vectors @var{x}, @var{y} of common length @math{N}, Kendall's
+## @tex
+## $\tau$
+## @end tex
+## @ifnottex
+## @var{tau}
+## @end ifnottex
+## is the correlation of the signs of all rank differences of
 ## @var{x} and @var{y}; i.e., if both @var{x} and @var{y} have distinct
 ## entries, then
 ##
 ## @tex
-## $$ \tau = {1 \over n(n-1)} \sum_{i,j} {\rm sign}(q_i-q_j) {\rm sign}(r_i-r_j) $$
+## $$ \tau = {1 \over N(N-1)} \sum_{i,j} {\rm sign}(q_i-q_j) \, {\rm sign}(r_i-r_j) $$
 ## @end tex
 ## @ifnottex
 ##
 ## @example
 ## @group
 ##          1
-## tau = -------   SUM sign (q(i) - q(j)) * sign (r(i) - r(j))
-##       n (n-1)   i,j
+## @var{tau} = -------   SUM sign (@var{q}(i) - @var{q}(j)) * sign (@var{r}(i) - @var{r}(j))
+##       N (N-1)   i,j
 ## @end group
 ## @end example
 ##
@@ -47,17 +59,24 @@
 ## $q_i$ and $r_i$
 ## @end tex
 ## @ifnottex
-## @var{q}(@var{i}) and @var{r}(@var{i})
+## @var{q}(i) and @var{r}(i)
 ## @end ifnottex
 ## are the ranks of @var{x} and @var{y}, respectively.
 ##
 ## If @var{x} and @var{y} are drawn from independent distributions,
-## Kendall's @var{tau} is asymptotically normal with mean 0 and variance
+## Kendall's
 ## @tex
-## ${2 (2n+5) \over 9n(n-1)}$.
+## $\tau$
 ## @end tex
 ## @ifnottex
-## @code{(2 * (2@var{n}+5)) / (9 * @var{n} * (@var{n}-1))}.
+## @var{tau}
+## @end ifnottex
+## is asymptotically normal with mean 0 and variance
+## @tex
+## ${2 (2N+5) \over 9N(N-1)}$.
+## @end tex
+## @ifnottex
+## @code{(2 * (2N+5)) / (9 * N * (N-1))}.
 ## @end ifnottex
 ##
 ## @code{kendall (@var{x})} is equivalent to @code{kendall (@var{x},

--- a/scripts/statistics/base/kurtosis.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/kurtosis.m	Wed Dec 27 23:38:25 2017 -0600
@@ -27,9 +27,9 @@
 ## @tex
 ## $$
 ## \kappa_1 = {{{1\over N}\,
-##          \sum_{i=1}^N (@var{x}_i - \bar{@var{x}})^4} \over \sigma^4},
+##          \sum_{i=1}^N (x_i - \bar{x})^4} \over \sigma^4},
 ## $$
-## where $N$ is the length of @var{x}, $\bar{@var{x}}$ its mean, and $\sigma$
+## where $N$ is the length of @var{x}, $\bar{x}$ its mean, and $\sigma$
 ## its (uncorrected) standard deviation.
 ## @end tex
 ## @ifnottex
@@ -65,6 +65,8 @@
 ## @end group
 ## @end example
 ##
+## where @math{N} is the length of the @var{x} vector.
+##
 ## @end ifnottex
 ## The bias-corrected kurtosis coefficient is obtained by replacing the sample
 ## second and fourth central moments by their unbiased versions.  It is an

--- a/scripts/statistics/base/mean.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/mean.m	Wed Dec 27 23:38:25 2017 -0600
@@ -27,13 +27,17 @@
 ##
 ## @tex
 ## $$ {\rm mean}(x) = \bar{x} = {1\over N} \sum_{i=1}^N x_i $$
+## where $N$ is the number of elements of @var{x}.
+##
 ## @end tex
 ## @ifnottex
 ##
 ## @example
-## mean (x) = SUM_i x(i) / N
+## mean (@var{x}) = SUM_i @var{x}(i) / N
 ## @end example
 ##
+## where @math{N} is the length of the @var{x} vector.
+##
 ## @end ifnottex
 ## If @var{x} is a matrix, compute the mean for each column and return them
 ## in a row vector.

--- a/scripts/statistics/base/meansq.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/meansq.m	Wed Dec 27 23:38:25 2017 -0600
@@ -27,18 +27,21 @@
 ## $$
 ## {\rm meansq} (x) = {\sum_{i=1}^N {x_i}^2 \over N}
 ## $$
-## where $\bar{x}$ is the mean value of $x$.
+## where $N$ is the number of elements of @var{x}.
+##
 ## @end tex
 ## @ifnottex
 ##
 ## @example
 ## @group
-## meansq (x) = 1/N SUM_i x(i)^2
+## meansq (@var{x}) = 1/N SUM_i @var{x}(i)^2
 ## @end group
 ## @end example
 ##
+## where @math{N} is the length of the @var{x} vector.
+##
 ## @end ifnottex
-## For matrix arguments, return a row vector containing the mean square
+## If @var{x} is a matrix, return a row vector containing the mean square
 ## of each column.
 ##
 ## If the optional argument @var{dim} is given, operate along this dimension.

--- a/scripts/statistics/base/median.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/median.m	Wed Dec 27 23:38:25 2017 -0600
@@ -22,27 +22,30 @@
 ## @deftypefnx {} {} median (@var{x}, @var{dim})
 ## Compute the median value of the elements of the vector @var{x}.
 ##
-## When the elements of @var{x} are sorted, the median is defined as
+## When the elements of @var{x} are sorted, say @code{@var{s} = sort (@var{x})},
+## the median is defined as
 ## @tex
 ## $$
 ## {\rm median} (x) =
-##   \cases{x(\lceil N/2\rceil), & $N$ odd;\cr
-##           (x(N/2)+x(N/2+1))/2, & $N$ even.}
+##   \cases{s(\lceil N/2\rceil), & $N$ odd;\cr
+##           (s(N/2)+s(N/2+1))/2, & $N$ even.}
 ## $$
+## where $N$ is the number of elements of @var{x}.
+##
 ## @end tex
 ## @ifnottex
 ##
 ## @example
 ## @group
-##               x(ceil(N/2))             N odd
-## median (x) =
-##              (x(N/2) + x((N/2)+1))/2   N even
+##              |  @var{s}(ceil(N/2))           N odd
+## median (@var{x}) = |
+##              | (@var{s}(N/2) + @var{s}(N/2+1))/2   N even
 ## @end group
 ## @end example
 ##
 ## @end ifnottex
 ## If @var{x} is of a discrete type such as integer or logical, then
-## the case of even N rounds up (or toward @code{true}).
+## the case of even @math{N} rounds up (or toward @code{true}).
 ##
 ## If @var{x} is a matrix, compute the median value for each column and
 ## return them in a row vector.

--- a/scripts/statistics/base/moment.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/moment.m	Wed Dec 27 23:38:25 2017 -0600
@@ -22,21 +22,26 @@
 ## @deftypefnx {} {} moment (@var{x}, @var{p}, @var{dim})
 ## @deftypefnx {} {} moment (@var{x}, @var{p}, @var{type}, @var{dim})
 ## @deftypefnx {} {} moment (@var{x}, @var{p}, @var{dim}, @var{type})
-## Compute the @var{p}-th central moment of the vector @var{x}.
+## Compute the @var{p}-th central moment of the vector @var{x}:
 ##
 ## @tex
 ## $$
 ## {\sum_{i=1}^N (x_i - \bar{x})^p \over N}
 ## $$
+## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of elements of @var{x}.
+##
+##
 ## @end tex
 ## @ifnottex
 ##
 ## @example
 ## @group
-## 1/N SUM_i (x(i) - mean(x))^p
+## 1/N SUM_i (@var{x}(i) - mean(@var{x}))^@var{p}
 ## @end group
 ## @end example
 ##
+## where @math{N} is the length of the @var{x} vector.
+##
 ## @end ifnottex
 ##
 ## If @var{x} is a matrix, return the row vector containing the @var{p}-th
@@ -64,7 +69,7 @@
 ##
 ## @example
 ## @group
-## 1/N SUM_i (abs (x(i) - mean(x)))^p
+## 1/N SUM_i (abs (@var{x}(i) - mean(@var{x})))^@var{p}
 ## @end group
 ## @end example
 ##
@@ -82,7 +87,7 @@
 ##
 ## @example
 ## @group
-## moment (x) = 1/N SUM_i x(i)^p
+## moment (@var{x}) = 1/N SUM_i @var{x}(i)^@var{p}
 ## @end group
 ## @end example
 ##
@@ -99,7 +104,7 @@
 ##
 ## @example
 ## @group
-## 1/N SUM_i ( abs (x(i)) )^p
+## 1/N SUM_i ( abs (@var{x}(i)) )^@var{p}
 ## @end group
 ## @end example
 ##

--- a/scripts/statistics/base/ols.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/ols.m	Wed Dec 27 23:38:25 2017 -0600
@@ -18,45 +18,61 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {} {[@var{beta}, @var{sigma}, @var{r}] =} ols (@var{y}, @var{x})
-## Ordinary least squares estimation.
+## Ordinary least squares (OLS) estimation.
 ##
 ## OLS applies to the multivariate model
 ## @tex
-## $y = x b + e$
-## with
-## $\bar{e} = 0$, and cov(vec($e$)) = kron ($s, I$)
+## $@var{y} = @var{x}\,@var{b} + @var{e}$
 ## @end tex
 ## @ifnottex
-## @w{@math{y = x*b + e}} with
-## @math{mean (e) = 0} and @math{cov (vec (e)) = kron (s, I)}.
+## @w{@math{@var{y} = @var{x}*@var{b} + @var{e}}}
 ## @end ifnottex
 ## where
 ## @tex
-## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix,
-## $b$ is a $k \times p$ matrix, and $e$ is a $t \times p$ matrix.
+## $@var{y}$ is a $t \times p$ matrix, $@var{x}$ is a $t \times k$ matrix,
+## $@var{b}$ is a $k \times p$ matrix, and $@var{e}$ is a $t \times p$ matrix.
 ## @end tex
 ## @ifnottex
-## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by @math{k}
-## matrix, @math{b} is a @math{k} by @math{p} matrix, and @math{e} is a
-## @math{t} by @math{p} matrix.
+## @math{@var{y}} is a @math{t}-by-@math{p} matrix, @math{@var{x}} is a
+## @math{t}-by-@math{k} matrix, @var{b} is a @math{k}-by-@math{p} matrix, and
+## @var{e} is a @math{t}-by-@math{p} matrix.
 ## @end ifnottex
 ##
-## Each row of @var{y} and @var{x} is an observation and each column a
-## variable.
+## Each row of @var{y} is a @math{p}-variate observation in which each column
+## represents a variable.  Likewise, the rows of @var{x} represent
+## @math{k}-variate observations or possibly designed values.  Furthermore,
+## the collection of observations @var{x} must be of adequate rank, @math{k},
+## otherwise @var{b} cannot be uniquely estimated.
+##
+## The observation errors, @var{e}, are assumed to originate from an
+## underlying @math{p}-variate distribution with zero mean and
+## @math{p}-by-@math{p} covariance matrix @var{S}, both constant conditioned
+## on @var{x}.  Furthermore, the matrix @var{S} is constant with respect to
+## each observation such that
+## @tex
+## $\bar{@var{e}} = 0$ and cov(vec(@var{e})) =  kron(@var{s},@var{I}).
+## @end tex
+## @ifnottex
+## @code{mean (@var{e}) = 0} and
+## @code{cov (vec (@var{e})) = kron (@var{s}, @var{I})}.
+## @end ifnottex
+## (For cases
+## that don't meet this criteria, such as autocorrelated errors, see
+## generalized least squares, gls, for more efficient estimations.)
 ##
 ## The return values @var{beta}, @var{sigma}, and @var{r} are defined as
 ## follows.
 ##
 ## @table @var
 ## @item beta
-## The OLS estimator for @math{b}.
+## The OLS estimator for matrix @var{b}.
 ## @tex
-## $beta$ is calculated directly via $(x^Tx)^{-1} x^T y$ if the matrix $x^Tx$ is
-## of full rank.
+## @var{beta} is calculated directly via $(@var{x}^T@var{x})^{-1} @var{x}^T
+## @var{y}$ if the matrix $@var{x}^T@var{x}$ 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.
+## @var{beta} is calculated directly via @code{inv (@var{x}'*@var{x}) * @var{x}' * @var{y}} if the
+## matrix @code{@var{x}'*@var{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}.
@@ -66,9 +82,7 @@
 ##
 ## @example
 ## @group
-## @var{sigma} = (@var{y}-@var{x}*@var{beta})'
-##   * (@var{y}-@var{x}*@var{beta})
-##   / (@var{t}-rank(@var{x}))
+## @var{sigma} = (@var{y}-@var{x}*@var{beta})' * (@var{y}-@var{x}*@var{beta}) / (@math{t}-rank(@var{x}))
 ## @end group
 ## @end example
 ##

--- a/scripts/statistics/base/prctile.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/prctile.m	Wed Dec 27 23:38:25 2017 -0600
@@ -24,7 +24,7 @@
 ## to the cumulative probability values, @var{p}, in percent.
 ##
 ## If @var{x} is a matrix, compute the percentiles for each column and return
-## them in a matrix, such that the i-th row of @var{y} contains the
+## them in a matrix, such that the i-th row of @var{q} contains the
 ## @var{p}(i)th percentiles of each column of @var{x}.
 ##
 ## If @var{p} is unspecified, return the quantiles for @code{[0 25 50 75 100]}.

--- a/scripts/statistics/base/quantile.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/quantile.m	Wed Dec 27 23:38:25 2017 -0600
@@ -37,7 +37,7 @@
 ##
 ## The methods available to calculate sample quantiles are the nine methods
 ## used by R (@url{http://www.r-project.org/}).  The default value is
-## @w{METHOD = 5}.
+## @w{@var{method} = 5}.
 ##
 ## Discontinuous sample quantile methods 1, 2, and 3
 ##
@@ -49,28 +49,72 @@
 ## @item Method 3: SAS definition: nearest even order statistic.
 ## @end enumerate
 ##
-## Continuous sample quantile methods 4 through 9, where p(k) is the linear
-## interpolation function respecting each methods' representative cdf.
+## Continuous sample quantile methods 4 through 9, where
+## @tex
+## $p(k)$
+## @end tex
+## @ifnottex
+## @var{p}(k)
+## @end ifnottex
+## is the linear
+## interpolation function respecting each method's representative cdf.
 ##
 ## @enumerate 4
-## @item Method 4: p(k) = k / n.  That is, linear interpolation of the
-## empirical cdf.
+## @item Method 4:
+## @tex
+## $p(k) = k / N$.
+## @end tex
+## @ifnottex
+## @var{p}(k) = k / N.
+## @end ifnottex
+## That is, linear interpolation of the empirical cdf, where @math{N} is the
+## length of @var{P}.
 ##
-## @item Method 5: p(k) = (k - 0.5) / n.  That is a piecewise linear function
-## where the knots are the values midway through the steps of the empirical
-## cdf.
+## @item Method 5:
+## @tex
+## $p(k) = (k - 0.5) / N$.
+## @end tex
+## @ifnottex
+## @var{p}(k) = (k - 0.5) / N.
+## @end ifnottex
+## That is, a piecewise linear function where the knots are the values midway
+## through the steps of the empirical cdf.
 ##
-## @item Method 6: p(k) = k / (n + 1).
+## @item Method 6:
+## @tex
+## $p(k) = k / (N + 1)$.
+## @end tex
+## @ifnottex
+## @var{p}(k) = k / (N + 1).
+## @end ifnottex
 ##
-## @item Method 7: p(k) = (k - 1) / (n - 1).
+## @item Method 7:
+## @tex
+## $p(k) = (k - 1) / (N - 1)$.
+## @end tex
+## @ifnottex
+## @var{p}(k) = (k - 1) / (N - 1).
+## @end ifnottex
 ##
-## @item Method 8: p(k) = (k - 1/3) / (n + 1/3).  The resulting quantile
-## estimates are approximately median-unbiased regardless of the distribution
-## of @var{x}.
+## @item Method 8:
+## @tex
+## $p(k) = (k - 1/3) / (N + 1/3)$.
+## @end tex
+## @ifnottex
+## @var{p}(k) = (k - 1/3) / (N + 1/3).
+## @end ifnottex
+## The resulting quantile estimates are approximately median-unbiased
+## regardless of the distribution of @var{x}.
 ##
-## @item Method 9: p(k) = (k - 3/8) / (n + 1/4).  The resulting quantile
-## estimates are approximately unbiased for the expected order statistics if
-## @var{x} is normally distributed.
+## @item Method 9:
+## @tex
+## $p(k) = (k - 3/8) / (N + 1/4)$.
+## @end tex
+## @ifnottex
+## @var{p}(k) = (k - 3/8) / (N + 1/4).
+## @end ifnottex
+## The resulting quantile estimates are approximately unbiased for the
+## expected order statistics if @var{x} is normally distributed.
 ## @end enumerate
 ##
 ## @nospell{Hyndman and Fan} (1996) recommend method 8.  Maxima, S, and R

--- a/scripts/statistics/base/skewness.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/skewness.m	Wed Dec 27 23:38:25 2017 -0600
@@ -27,9 +27,9 @@
 ## @tex
 ## $$
 ## {\rm skewness} (@var{x}) = {{{1\over N}\,
-##          \sum_{i=1}^N (@var{x}_i - \bar{@var{x}})^3} \over \sigma^3},
+##          \sum_{i=1}^N (x_i - \bar{x})^3} \over \sigma^3},
 ## $$
-## where $N$ is the length of @var{x}, $\bar{@var{x}}$ its mean and $\sigma$
+## where $N$ is the length of @var{x}, $\bar{x}$ its mean and $\sigma$
 ## its (uncorrected) standard deviation.
 ## @end tex
 ## @ifnottex
@@ -52,7 +52,7 @@
 ## @tex
 ## $$
 ## {\rm skewness} (@var{x}) = {\sqrt{N (N - 1)} \over N - 2} \times \,
-##   {{{1 \over N} \sum_{i=1}^N (@var{x}_i - \bar{@var{x}})^3} \over \sigma^3}
+##   {{{1 \over N} \sum_{i=1}^N (x_i - \bar{x})^3} \over \sigma^3}
 ## $$
 ## @end tex
 ## @ifnottex
@@ -65,6 +65,8 @@
 ## @end group
 ## @end example
 ##
+## where @math{N} is the length of the @var{x} vector.
+##
 ## @end ifnottex
 ## The adjusted skewness coefficient is obtained by replacing the sample second
 ## and third central moments by their bias-corrected versions.

--- a/scripts/statistics/base/spearman.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/spearman.m	Wed Dec 27 23:38:25 2017 -0600
@@ -20,14 +20,39 @@
 ## @deftypefn  {} {} spearman (@var{x})
 ## @deftypefnx {} {} spearman (@var{x}, @var{y})
 ## @cindex Spearman's Rho
-## Compute Spearman's rank correlation coefficient @var{rho}.
+## Compute Spearman's rank correlation coefficient
+## @tex
+## $\rho$.
+## @end tex
+## @ifnottex
+## @var{rho}.
+## @end ifnottex
+##
+## For two data vectors @var{x} and @var{y}, Spearman's
+## @tex
+## $\rho$
+## @end tex
+## @ifnottex
+## @var{rho}
+## @end ifnottex
+## is the correlation coefficient of the ranks of @var{x} and @var{y}.
 ##
-## For two data vectors @var{x} and @var{y}, Spearman's @var{rho} is the
-## correlation coefficient of the ranks of @var{x} and @var{y}.
-##
-## If @var{x} and @var{y} are drawn from independent distributions, @var{rho}
-## has zero mean and variance @code{1 / (n - 1)}, and is asymptotically
-## normally distributed.
+## If @var{x} and @var{y} are drawn from independent distributions,
+## @tex
+## $\rho$
+## @end tex
+## @ifnottex
+## @var{rho}
+## @end ifnottex
+## has zero mean and variance
+## @tex
+## $1 / (N - 1)$,
+## @end tex
+## @ifnottex
+## @code{1 / (N - 1)},
+## @end ifnottex
+## where @math{N} is the length of the @var{x} and @var{y} vectors, and is
+## asymptotically normally distributed.
 ##
 ## @code{spearman (@var{x})} is equivalent to
 ## @code{spearman (@var{x}, @var{x})}.

--- a/scripts/statistics/base/std.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/std.m	Wed Dec 27 23:38:25 2017 -0600
@@ -27,18 +27,18 @@
 ## $$
 ## {\rm std} (x) = \sigma = \sqrt{{\sum_{i=1}^N (x_i - \bar{x})^2 \over N - 1}}
 ## $$
-## where $\bar{x}$ is the mean value of $x$ and $N$ is the number of elements.
+## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of elements of @var{x}.
 ## @end tex
 ## @ifnottex
 ##
 ## @example
 ## @group
-## std (x) = sqrt ( 1/(N-1) SUM_i (x(i) - mean(x))^2 )
+## std (@var{x}) = sqrt ( 1/(N-1) SUM_i (@var{x}(i) - mean(@var{x}))^2 )
 ## @end group
 ## @end example
 ##
 ## @noindent
-## where @math{N} is the number of elements.
+## where @math{N} is the number of elements of the @var{x} vector.
 ## @end ifnottex
 ##
 ## If @var{x} is a matrix, compute the standard deviation for each column and

--- a/scripts/statistics/base/var.m	Wed Jan 03 21:49:38 2018 -0800
+++ b/scripts/statistics/base/var.m	Wed Dec 27 23:38:25 2017 -0600
@@ -27,16 +27,20 @@
 ## $$
 ## {\rm var} (x) = \sigma^2 = {\sum_{i=1}^N (x_i - \bar{x})^2 \over N - 1}
 ## $$
-## where $\bar{x}$ is the mean value of $x$.
+## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of
+## elements of @var{x}.
+##
 ## @end tex
 ## @ifnottex
 ##
 ## @example
 ## @group
-## var (x) = 1/(N-1) SUM_i (x(i) - mean(x))^2
+## var (@var{x}) = 1/(N-1) SUM_i (@var{x}(i) - mean(@var{x}))^2
 ## @end group
 ## @end example
 ##
+## where @math{N} is the length of the @var{x} vector.
+##
 ## @end ifnottex
 ## If @var{x} is a matrix, compute the variance for each column and return
 ## them in a row vector.
@@ -53,8 +57,8 @@
 ##   normalizes with @math{N}, this provides the second moment around the mean
 ## @end table
 ##
-## If @math{N==1} the value of @var{opt} is ignored and normalization by
-## @math{N} is used.
+## If @math{N} is equal to 1 the value of @var{opt} is ignored and
+## normalization by @math{N} is used.
 ##
 ## If the optional argument @var{dim} is given, operate along this dimension.
 ## @seealso{cov, std, skewness, kurtosis, moment}