(* :Name: Statistics`MultiDescriptiveStatistics` *)

(* :Title: Multivariate Descriptive Statistics *)

(* :Author: E. C. Martin *)

(* :Summary:
This package computes descriptive statistics (location, dispersion, shape,
and association statistics) for a sample represented as a data matrix.
The data matrix is a list of independent identically distributed
vector-valued or multivariate observations.
*)

(* :Context: Statistics`MultiDescriptiveStatistics` *)

(* :Package Version: 1.0 *)

(* :Copyright: Copyright 1993-1999  by Wolfram Research, Inc. *)

(* :History:
Original version by ECM (Wolfram Research), September 1993.
*)

(* :Source:
	Multivariate Analysis, Mardia, Kent, Bibby, Academic Press, 1979.
	A Survey of Multidimensional Medians, C. G. Small, International
		Statistical Review, 58, 3, pp. 263-277, 1990.
	Robust Statistics, P. J. Huber, John Wiley & Sons, Inc., 1980.
	Robust Estimates, Residuals, and Outlier Detection with
		Multiresponse Data, R. Gnanadesikan & J. R. Kettenring,
		Biometrics 28, March 1972.
	Rank Correlation Methods, M. G. Kendall, Hafner Publishing Co., 1955.
*)

(* :Warnings: 
	Functions based on ConvexHull currently only available for 
	dimension 2.
	Rules attached to Graphics and Graphics3D. *)
(* :Limitations: 
	Functions based on ConvexHull currently only available for 
	dimension 2. *)

(* :Mathematica Version: 3.0 *)



BeginPackage["Statistics`MultiDescriptiveStatistics`",
	(* needed for Mean, Quantile, etc: *)
	"Statistics`DescriptiveStatistics`", 
	(* needed for ChiSquareDistribution: *)
	"Statistics`NormalDistribution`",
	(* needed for CovarianceMatrix, CorrelationMatrix, and Ellipsoid: *)
	"Statistics`Common`MultivariateCommon`",
	(* needed for ConvexHull: *)
	"DiscreteMath`ComputationalGeometry`"]

Unprotect[Mean, GeometricMean, HarmonicMean, RootMeanSquare,
	TrimmedMean, Quantile, InterpolatedQuantile, Median, Mode, Quartiles,
	SampleRange, Variance, VarianceMLE, VarianceOfSampleMean,
	StandardDeviation, StandardDeviationMLE, StandardErrorOfSampleMean,
 	MeanDeviation, MedianDeviation, InterquartileRange,
	QuartileDeviation, CentralMoment,
	Skewness, PearsonSkewness1, PearsonSkewness2,
	QuartileSkewness, Kurtosis, KurtosisExcess];

(* Add multivariate usage to Mode.  Mode is treated separately
	since it is a system synbol. *)
If[StringPosition[Mode::usage, "matrix"] === {},

Mode::usage = Mode::usage <> " " <>
"Mode[matrix] gives the mode(s) of each column in matrix.";
]


(* Add multivariate usage to various univariate functions. *)
If[StringPosition[Mean::usage, "matrix"] === {},

ExpectedValue::usage = ExpectedValue::usage <> " " <>
"ExpectedValue[f, matrix] gives the expected value of the pure function f
of p variables with respect to the p-dimensional sample distribution of the
n observations in the n by p matrix.
ExpectedValue[f, matrix, {x1, x2, ..., xp}] gives the expected value of the
function f of x1, x2, ..., xp with respect to the sample distribution of
matrix.";

Mean::usage = Mean::usage <> " " <>
"Mean[matrix] gives the mean of each column in matrix.";

GeometricMean::usage = GeometricMean::usage <> " " <>
"GeometricMean[matrix] gives the geometric mean of each column in matrix.";

HarmonicMean::usage = HarmonicMean::usage <> " " <>
"HarmonicMean[matrix] gives the geometric mean of each column in matrix.";

RootMeanSquare::usage = RootMeanSquare::usage <> " " <>
"RootMeanSquare[matrix] gives the root mean square of each column in matrix.";

TrimmedMean::usage = TrimmedMean::usage <> " " <>
"TrimmedMean[matrix, f] gives the mean of each column in matrix, with a fraction
f of entries at each end dropped.  TrimmedMean[matrix, {f1, f2}] finds the mean
with fractions f1 and f2 trimmed off.";

Quantile::usage = Quantile::usage <> " " <>
"Quantile[matrix, q] gives the q-th quantile of each column in matrix.";

InterpolatedQuantile::usage = InterpolatedQuantile::usage <> " " <>
"InterpolatedQuantile[matrix, q] gives the q-th interpolated quantile of each
column in matrix.";

Median::usage = Median::usage <> " " <>
"Median[matrix] gives the median of each column in matrix (the coordinate-wise
median).";

Quartiles::usage = Quartiles::usage <> " " <>
"Quartiles[matrix] gives the quartiles of each column in matrix.";

SampleRange::usage = SampleRange::usage <> " " <>
"SampleRange[matrix] gives the range of each column in matrix.";

Variance::usage = Variance::usage <> " " <>
"Variance[matrix] gives an unbiased estimate of the variance of each column in
matrix.";

VarianceMLE::usage = VarianceMLE::usage <> " " <>
"VarianceMLE[matrix] gives the maximum likelihood estimate of the variance of 
each column in matrix.";

VarianceOfSampleMean::usage = VarianceOfSampleMean::usage <> " " <>
"VarianceOfSampleMean[matrix] returns an unbiased estimate of the variance
of the sample mean of each column in matrix.";

StandardDeviation::usage = StandardDeviation::usage <> " " <>
"StandardDeviation[matrix] gives an unbiased estimate of the standard deviation
of each column in matrix.";

StandardDeviationMLE::usage = StandardDeviationMLE::usage <> " " <>
"StandardDeviationMLE[matrix] returns the maximum likelihood estimate of the
standard deviation of each column in matrix.";

StandardErrorOfSampleMean::usage = StandardErrorOfSampleMean::usage <> " " <>
"StandardErrorOfSampleMean[matrix] returns an unbiased estimate of the
standard error of the sample mean for each column in matrix.";

MeanDeviation::usage = MeanDeviation::usage <> " " <>
"MeanDeviation[matrix] gives the mean absolute deviation about the mean for
each column in matrix.";

MedianDeviation::usage = MedianDeviation::usage <> " " <>
"MedianDeviation[matrix] gives the median absolute deviation about the median
for each column in matrix.";

InterquartileRange::usage = InterquartileRange::usage <> " " <>
"InterquartileRange[matrix] gives the difference between the upper and lower
quartiles for each column in matrix.";

QuartileDeviation::usage = QuartileDeviation::usage <> " " <>
"QuartileDeviation[matrix] gives the quartile deviation for each column in
matrix.";

CentralMoment::usage = CentralMoment::usage <> " " <>
"CentralMoment[matrix, r] gives the r-th central moment for each column in
matrix.  CentralMoment[matrix, {r1, ..., rp}] gives an r-th central moment
of the n p-dimensional vectors in the n x p matrix, where r = r1 + ... + rp.";

Skewness::usage = Skewness::usage <> " " <>
"Skewness[matrix] gives the coefficient of skewness for each column in matrix.";

PearsonSkewness1::usage = PearsonSkewness1::usage <> " " <>
"PearsonSkewness1[matrix] gives Pearson's first coefficient of skewness for each
column in matrix.";

PearsonSkewness2::usage = PearsonSkewness2::usage <> " " <>
"PearsonSkewness2[matrix] gives Pearson's second coefficient of skewness for
each column in matrix.";

QuartileSkewness::usage = QuartileSkewness::usage <> " " <>
"QuartileSkewness[matrix] gives the quartile coefficient of skewness for each
column in matrix.";

Kurtosis::usage = Kurtosis::usage <> " " <>
"Kurtosis[matrix] gives the kurtosis coefficient for each column in matrix.";

KurtosisExcess::usage = KurtosisExcess::usage <> " " <>
"KurtosisExcess[matrix] gives the  kurtosis excess for each column in matrix.";

ZeroMean::usage = ZeroMean::usage <> " " <>
"ZeroMean[matrix] subtracts the mean from each column of matrix.";

Standardize::usage = Standardize::usage <> " " <>
"Standardize[matrix] standardizes each column of matrix to have zero mean and
unit variance.  Standardize[matrix, Decorrelate->True] applies the
Mahalanobis transformation so that the resulting p-variate data have covariance
equal to the p-dimensional identity matrix."

] (* end If *)

Decorrelate::usage =
"Decorrelate is an option of Standardize indicating whether the data should
be transformed so that the covariances are zero."

(* multivariate location statistics *)

SpatialMedian::usage =
"SpatialMedian[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] estimates the
p-dimensional median to be that vector minimizing the sum of the Euclidean
distances between the vector and the n p-dimensional points in the sample."

SimplexMedian::usage =
"SimplexMedian[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] estimates the
p-dimensional median to be that vector minimizing the sum of the volumes
of the p-dimensional simplices the vector forms with all possible combinations
of p members from the sample."

ConvexHullMedian::usage =	
"ConvexHullMedian[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] estimates the
p-dimensional median to be the mean of the p-dimensional vectors lying on the
innermost layer of the convex layers of the n p-dimensional points."

EstimateDOF::usage =
"EstimateDOF is an option to several multivariate
descriptive statistic functions that give robust results by trimming
outliers.  EstimateDOF -> False (default) means that only the
estimate itself is returned.  EstimateDOF -> True means that
both the estimate and the degrees of freedom in the estimate (number
of points contributing to the estimate) are returned in the form
{estimate, EstimateDOF -> dof}."

MultivariateTrimmedMean::usage = 
"MultivariateTrimmedMean[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}, f] gives the
mean of the p-dimensional vectors, with a fraction f of the most outlying
vectors dropped.  When f = 0, MultivariateTrimmedMean gives the mean.
As f -> 1.0, MultivariateTrimmedMean approaches the multivariate median
ConvexHullMedian."

MultivariateMode::usage = 
"MultivariateMode[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the mode of
the n p-dimensional vectors.  A list of modes is returned if the data is
bimodal, trimodal, or multimodal."

EllipsoidQuantile::usage =
"EllipsoidQuantile[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}, q] gives the locus
of the qth quantile of the p-variate data, where the data have been ordered
using ellipsoids centered on Mean[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}].
The fraction of the data lying inside of this locus is q."

EllipsoidQuartiles::usage =
"EllipsoidQuartiles[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives a list of 
the locii of the quartiles (q = .25, .50, .75) of the p-variate data, where 
the data have been ordered using ellipsoids centered on 
Mean[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}]."

PolytopeQuantile::usage =
"PolytopeQuantile[{{x11, x12}, ..., {xn1, xn2}}, q] gives the locus of the
qth quantile of the bivariate data, where the data have been ordered using
convex hulls centered on ConvexHullMedian[{{x11, x12}, ..., {xn1, xn2}}].
The fraction of the data lying inside of this locus is q."

PolytopeQuartiles::usage =
"PolytopeQuartiles[{{x11, x12}, ..., {xn1, xn2}}] gives a list of the locii
of the quartiles (q = .25, .50, .75) of the bivariate data, where the data
have been ordered using convex hulls centered on 
ConvexHullMedian[{{x11, x12}, ..., {xn1, xn2}}]."

(* bivariate dispersion statistics *)

Covariance::usage =
"Covariance[{x1, ...., xn}, {y1, ...., yn}] gives the covariance between the
x and y variables.  Division by n-1 (rather than n) is used, giving an
unbiased estimate of the population covariance.  Covariance[xlist, ylist,
ScaleMethod -> method] computes covariance using a measure of scale other than
StandardDeviation (i.e., MeanDeviation, MedianDeviation, QuartileDeviation);
scale is standardized to give unity when the data is normal with zero mean
and unit variance."

CovarianceMLE::usage =
"CovarianceMLE[{x1, ...., xn}, {y1, ...., yn}] gives the covariance between 
the x and y variables.  Division by n (rather than n-1) is used, giving a
maximum likelihood estimate of the population covariance."

ScaleMethod::usage =
"ScaleMethod is an option of multivariate descriptive statistic functions
indicating which estimate of scale (MeanDeviation, MedianDeviation,
QuartileDeviation) is used.  The default is ScaleMethod->StandardDeviation."

(* matrix-valued multivariate dispersion statistics *)

CovarianceMatrixMLE::usage =
"CovarianceMatrixMLE[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
p x p covariance matrix of the n p-dimensional vectors.  Division by n
(rather than n-1) is used, giving a maximum likelihood estimate of the
population covariance (use CovarianceMatrix for an unbiased estimate).
CovarianceMatrixMLE[{{x11, ..., x1p}, ..., {xn1, ..., xnp}},
{{y11, ..., y1q}, ..., {yn1, ..., ynq}}] gives the p x q covariance matrix 
between the n p-dimensional vectors and the n q-dimensional vectors."

CovarianceMatrixOfSampleMean::usage =
"CovarianceMatrixOfSampleMean[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] returns 
an unbiased estimate of the p x p covariance matrix of the p-dimensional sample
mean vector, using the vectors {xi1, ..., xip}, i = 1, ..., n, as a sample from 
the population."

DispersionMatrix::usage =
"DispersionMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] finds a p x p  
dispersion matrix describing the n p-dimensional vectors, by computing
a covariance matrix from only those points lying inside the convex hull of
the points."

(* bivariate association statistics *)

Correlation::usage =
"Correlation[{x1, ...., xn}, {y1, ...., yn}] gives the linear correlation 
coefficient between the x and y variables.  Correlation[xlist, ylist,
ScaleMethod -> method] computes correlation using a measure of scale other than
StandardDeviation (i.e., MeanDeviation, MedianDeviation, QuartileDeviation)."

SpearmanRankCorrelation::usage =
"SpearmanRankCorrelation[{x1, ...., xn}, {y1, ...., yn}] gives Spearman's
rank correlation coefficient (termed rho-b) between the x and y variables.
(The alternatives rho-a and rho-b are different approaches for handling
ties in rankings.  The absolute value of rho-b is greater than the 
absolute value of rho-a in the case of ties.)"

KendallRankCorrelation::usage =
"KendallRankCorrelation[{x1, ...., xn}, {y1, ...., yn}] gives Kendall's
rank correlation coefficient (termed tau-b) between the x and y variables.
(The alternatives tau-a and tau-b are different approaches for handling
ties in rankings.  The absolute value of tau-b is greater than the 
absolute value of tau-a in the case of ties.)"

(* multivariate association statistics *)

AssociationMatrix::usage =
"AssociationMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] finds a p x p  
association matrix describing the n p-dimensional vectors, by computing
a correlation matrix from only those points lying inside the convex hull of
the points."

(* scalar-valued multivariate dispersion statistics *)

GeneralizedVariance::usage =
"GeneralizedVariance[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
generalized variance of the n p-dimensional vectors.  This is equivalent to
the determinant of the covariance matrix, or the product of the variances of the
principal components of the data."

TotalVariation::usage =
"TotalVariation[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the total
variation of the n p-dimensional vectors.  This is equivalent to the trace of
the covariance matrix, or the sum of the variances of the principal components
of the data."

ConvexHullArea::usage =
"ConvexHullArea[{{x1, y1}, ..., {xn, yn}}] gives the area
of the convex hull of the bivariate data."
(* This should be extended to dimensions > 2 as soon as convex hull for
dimensions > 2 is implemented. *)

MultivariateMeanDeviation::usage =
"MultivariateMeanDeviation[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
scalar mean of the Euclidean distances between the p-variate mean and the
p-variate data."

MultivariateMedianDeviation::usage =
"MultivariateMedianDeviation[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
scalar median of the Euclidean distances between the p-variate median and the
p-variate data."

MedianMethod::usage =
"MedianMethod is an option of multivariate descriptive statistic
functions indicating which estimate of multivariate median (Median,
SpatialMedian, SimplexMedian, or ConvexHullMedian) is used."

(* multivariate shape statistics *)

MultivariateSkewness::usage =
"MultivariateSkewness[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
multivariate coefficient of skewness for the p-variate data.  A value close
to zero indicates elliptical symmetry."

MultivariatePearsonSkewness1::usage =
"MultivariatePearsonSkewness1[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives
9 (mean-mode)'.sigma^(-1).(mean-mode), where mean is the mean vector, mode is
the mode vector, and sigma is the covariance matrix.  This is a multivariate
extension of PearsonSkewness1 which gives 3(mean-mode)/sigma, where sigma is
the standard deviation."

MultivariatePearsonSkewness2::usage =
"MultivariatePearsonSkewness1[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives
9 (mean-median)'.sigma^(-1).(mean-median), where mean is the mean vector, mode 
is the mode vector, and sigma is the covariance matrix.  This is a multivariate
extension of PearsonSkewness2 which gives 3(mean-median)/sigma, where sigma is
the standard deviation."

MultivariateKurtosis::usage =
"MultivariateKurtosis[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
multivariate kurtosis coefficient for the p-variate data.  A value close to
p*(p+2) indicates multinormality."

MultivariateKurtosisExcess::usage =
"MultivariateKurtosisExcess[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
multivariate kurtosis excess for the p-variate data.  A value close to zero
indicates multinormality."

(* multivariate data transformations *)

PrincipalComponents::usage =
"PrincipalComponents[matrix] transforms each row of p variables to p principal
components.  This is a set of uncorrelated variables, arranged in order of
decreasing variance."


(* updating usage of shared symbols *)

If[StringPosition[CovarianceMatrix::usage, "n p-dimensional"] === {},

(* CovarianceMatrix::usage exists, but multivariate usage does not. *)
CovarianceMatrix::usage = CovarianceMatrix::usage <> " " <>
"CovarianceMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
p x p covariance matrix of the n p-dimensional vectors.  Division by n-1
(rather than n) is used, giving an unbiased estimate of the population
covariance (use CovarianceMatrixMLE for a maximum likelihood estimate).
CovarianceMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}},
{{y11, ..., y1q}, ..., {yn1, ..., ynq}}] gives the p x q covariance matrix 
between the n p-dimensional vectors and the n q-dimensional vectors.";

CorrelationMatrix::usage = CorrelationMatrix::usage <> " " <>
"CorrelationMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
p x p correlation matrix of the n p-dimensional vectors.
CorrelationMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}},
{{y11, ..., y1q}, ..., {yn1, ..., ynq}}] gives the p x q correlation matrix
between the n p-dimensional vectors and the n q-dimensional vectors.",

(* CovarianceMatrix::usage does not exist. *)
CovarianceMatrix::usage =
"CovarianceMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
p x p covariance matrix of the n p-dimensional vectors.  Division by n-1
(rather than n) is used, giving an unbiased estimate of the population
covariance (use CovarianceMatrixMLE for a maximum likelihood estimate).
CovarianceMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}},
{{y11, ..., y1q}, ..., {yn1, ..., ynq}}] gives the p x q covariance matrix 
between the n p-dimensional vectors and the n q-dimensional vectors.";

CorrelationMatrix::usage =
"CorrelationMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}] gives the
p x p correlation matrix of the n p-dimensional vectors.
CorrelationMatrix[{{x11, ..., x1p}, ..., {xn1, ..., xnp}},
{{y11, ..., y1q}, ..., {yn1, ..., ynq}}] gives the p x q correlation matrix
between the n p-dimensional vectors and the n q-dimensional vectors."

]


Begin["`Private`"]

(* ========== multivariate extensions of univariate functions =========== *)

makemultivariate1[f_] :=
  (
  f[m_?MatrixQ, args___] := Module[{v}, v /;
      FreeQ[v = Map[f[#, args]&, Transpose[m]], f]]
  )

Scan[makemultivariate1, 
 {Mean, GeometricMean, HarmonicMean, RootMeanSquare, TrimmedMean, Quantile,
  InterpolatedQuantile, Median,
  Quartiles, SampleRange, Variance, VarianceMLE, VarianceOfSampleMean,
  StandardDeviation, StandardDeviationMLE, StandardErrorOfSampleMean,
  MeanDeviation, MedianDeviation, InterquartileRange, QuartileDeviation,
  CentralMoment, Skewness, PearsonSkewness1, PearsonSkewness2,
  QuartileSkewness, Kurtosis, KurtosisExcess}
]

ExpectedValue[f_Function, matrix_?MatrixQ] :=
	Apply[Plus, Map[Apply[f, #]&, matrix]] / Length[matrix] 

ExpectedValue[f_, matrix_?MatrixQ, x_List] := Apply[Plus,
   Map[Apply[Function[Evaluate[x], Evaluate[f]], #]&, matrix]] / Length[matrix]

ZeroMean[m_?MatrixQ] :=
  Module[{m1},
     Transpose[m1] /; FreeQ[m1 = Map[ZeroMean, Transpose[m]], ZeroMean]
  ]

Options[Standardize] = {Decorrelate -> False, MLE -> False}

Standardize[m_?MatrixQ, opts___] :=
  Module[{m1, mle, s, sinv, sqrtsinv, mean},
     m1 /; 
     (mle = MLE /. {opts} /. Options[Standardize];
      If[TrueQ[Decorrelate /. {opts} /. Options[Standardize]],
	 (* decorrelate *)
	 If[If[TrueQ[!mle],
	       FreeQ[s = CovarianceMatrix[m], CovarianceMatrix],
	       FreeQ[s = CovarianceMatrixMLE[m], CovarianceMatrixMLE]
            ] &&
	    FreeQ[sinv = Inverse[s], Inverse] &&
	    FreeQ[svd = SingularValues[sinv], SingularValues] &&
	    ({u, w, v} = svd;
	     p = Length[m[[1]]];
	     Dimensions[u] == {p, p} && u==v),

	      sqrtsinv = DiagonalMatrix[Sqrt[w]].u;
	      mean = Mean[m];
	      m1 = Map[(sqrtsinv.(#-mean))&, m],
	      m1 = $Failed
         ],
	 (* do not decorrelate *)
	 If[FreeQ[m1 = Map[Standardize[#, MLE->mle]&, Transpose[m]],
		  Standardize],	
	    m1 = Transpose[m1],
	    m1 = $Failed
	 ]
      ];
      m1 =!= $Failed
     )
  ]


(* ============================= MultivariateMode ========================== *)
(* ==================  multivariate estimate of location =================== *)

MultivariateMode[m_?MatrixQ] :=
  Module[{u, c, maxc, ms},
    u = Union[m];
    c = Map[Count[m, #]&, u];
    maxc = Max[c];
    ms = Cases[Transpose[{u, c}], {val_, maxc}->val];
    If[maxc > 1,
	If[Length[ms] == 1, ms[[1]], ms],
	{}]
  ]  


(* =========================== EllipsoidQuantile =========================== *)
(* ==================  multivariate estimate of location =================== *)

Options[EllipsoidQuantile] = {MLE -> False}

EllipsoidQuantile[m_?MatrixQ, q_, opt___] :=
  Module[{quantile, s, sinv, mean, dist, layer, n, qlevel},
    (
    quantile
    ) /; If[TrueQ[MLE /. {opt} /. Options[EllipsoidQuantile]], 
	    FreeQ[s = CovarianceMatrixMLE[m], CovarianceMatrixMLE],
	    FreeQ[s = CovarianceMatrix[m], CovarianceMatrix]] &&
         FreeQ[sinv = Inverse[s], Inverse] &&
         (mean = Mean[m];
	  {dist, layer} = ellipticLayer[Map[Sqrt[(#-mean).sinv.(#-mean)]&, m]]; 
          (* qlevel = {1., (n-n1)/n, (n-n1-n2)/n, ...} corresponding to the 
             fact that the 1st layer contains 100% of the points, the 2nd 
             contains (n-n1)/n of the points, etc. *)
	  n = Length[m];
          qlevel = Drop[N[(n - FoldList[Plus, 0, Map[Length, layer]]) / n], -1];
          quantile = iEllipsoidQuantile[s, mean, qlevel, dist, q];
	  quantile =!= $Failed
         )
  ] /; Head[q]=!=List && NumberQ[N[q]] && FreeQ[N[q], Complex] && 0<=N[q]<=1

iEllipsoidQuantile[s_, mu_, qlevel_, dist_, q_] :=
  Module[{outer, inner, do, di, qi, distance},
    (* locate quantile between layers for a local fit *)
    outer = 1;
    Scan[If[q > #, Return[outer], outer++]&, Drop[qlevel, 1]];
    do = dist[[outer]];
    If[outer == Length[dist],
	{qi, di} = {0, 0},
	inner = outer+1;
	qi = qlevel[[inner]]; 
	di = dist[[inner]]
    ];
    distance = di + (q-qi)/(qlevel[[outer]]-qi) (do-di);
    ellipsoidalLocus[mu, distance^2 s]
  ]

(* Case where q is a list of quantiles. *)
EllipsoidQuantile[m_?MatrixQ, q_?VectorQ, opt___] :=
  Module[{quantiles, quantile, s, sinv, mean, dist, layer, n, qlevel, scan},
    (
    quantiles
    ) /; If[TrueQ[MLE /. {opt} /. Options[EllipsoidQuantile]],
	    FreeQ[s = CovarianceMatrixMLE[m], CovarianceMatrixMLE],
	    FreeQ[s = CovarianceMatrix[m], CovarianceMatrix]] &&
         FreeQ[sinv = Inverse[s], Inverse] &&
         (mean = Mean[m];
	  {dist, layer} = ellipticLayer[Map[Sqrt[(#-mean).sinv.(#-mean)]&, m]]; 
          (* qlevel = {1., (n-n1)/n, (n-n1-n2)/n, ...} corresponding to the 
             fact that the 1st layer contains 100% of the points, the 2nd 
             contains (n-n1)/n of the points, etc. *)
	  n = Length[m];
          qlevel = Drop[N[(n - FoldList[Plus, 0, Map[Length, layer]]) / n], -1];
	  quantiles = {};
	  scan = Scan[(quantile =
			 iEllipsoidQuantile[s, mean, qlevel, dist, #];
		       If[quantile =!= $Failed,
			  quantiles = Append[quantiles, quantile],
			  Return[$Failed]])&,
		      q];
          scan =!= $Failed)
  ] /; Apply[And,
	 Map[(NumberQ[#] && FreeQ[#, Complex] && 0<=#<=1)&, N[q]]]

(* provides list of distances together with indices of points that are
	that distance from the mean; distances listed from largest to
	smallest; {{d1, ..., dn}, {ilist1, ..., ilist2}} *)	
ellipticLayer[d_] :=
   Module[{distance, distinct},
     distance = Transpose[{d, Range[Length[d]]}];
     distinct = Sort[Union[Map[#[[1]]&, distance]], Greater];
     distinct = Map[Cases[distance, {#, _}]&, distinct];
     Transpose[ Map[Function[{y}, {y[[1, 1]], Map[#[[2]]&, y]}],
                    distinct] ]
   ]

(* used by EllipsoidQuantile *)
ellipsoidalLocus[mu_, sigma_] :=
  Module[{esystem, esystemT, sorted, r, dir},
    (* (x-mu).Inverse[sigma].(x-mu) == 1 *)
    If[!FreeQ[esystem = Eigensystem[sigma], Eigensystem], Return[$Failed]];
    (* radii are square roots of eigenvalues *)
    esystemT = Transpose[MapAt[Sqrt[Chop[#]]&, esystem, 1]];
    (* sort semi-axes from largest to smallest *)
    sorted = Sort[esystemT, #1[[1]] > #2[[1]]&];
    {r, dir} = Transpose[sorted];
    Ellipsoid[mu, r, dir]
  ]


(* =========================== EllipsoidQuartiles ========================== *)
(* ==================  multivariate estimate of location =================== *)

Options[EllipsoidQuartiles] = {MLE -> False}

EllipsoidQuartiles[m_?MatrixQ, opt___] :=
  EllipsoidQuantile[m, {.25, .50, .75}, opt] 


(* ========= SpatialMedian multivariate robust measure of location ========= *)
Options[SpatialMedian] =
	{AccuracyGoal -> Automatic, 
	 MaxIterations -> 60, PrecisionGoal -> Automatic,
	 WorkingPrecision -> Precision[1.]}


SpatialMedian::invpar =
"Warning: `` `` must be a positive integer.  Default value assumed."

SpatialMedian::findmin =
"Minimization of objective function using FindMinimum failed."

SpatialMedian[m_?MatrixQ, opts___] :=
   Module[{median = iSpatialMedian[m, opts]},
	median /; median =!= $Failed
   ]

iSpatialMedian[m_, opts___] :=
   Module[{n, p, coordvar, coord, 
	   ag, maxiter, pg, wp, 
	   marginalMedian, L1marg, mean, L1mean, init, L1init,
	   result, distlist, args, 
	   FMresult, L1result},
	{n, p} = Dimensions[m];
	coordvar = Table[Unique[coord], {p}]; 
	{ag, maxiter, pg, wp} =
	    {AccuracyGoal, MaxIterations, PrecisionGoal, WorkingPrecision} /.
		{opts} /. Options[SimplexMedian];
	If[!IntegerQ[maxiter] || maxiter < 1,
	   Message[SpatialMedian::invpar, MaxIterations, maxiter];
	   maxiter = MaxIterations /. Options[SpatialMedian]
        ];
	If[!IntegerQ[wp] || wp < 1,
	   Message[SpatialMedian::invpar, WorkingPrecision, wp];
	   wp = WorkingPrecision /. Options[SpatialMedian]
        ];

	marginalMedian = Median[m];
	L1marg = L1[m, marginalMedian];
	mean = Mean[m];
	L1mean = L1[m, mean];
	{init, L1init} = If[L1marg < L1mean,
				{marginalMedian, L1marg},
				{mean, L1mean}];

	distlist = Map[Sqrt[Apply[Plus, (coordvar - #)^2]]&, m];

	distlist = Map[If[Precision[#] < wp,
			   SetPrecision[#, wp],
			   #]&,
		        distlist];
	args = Join[ { Apply[Plus, distlist] },
		       Transpose[{coordvar, init}],
		     { AccuracyGoal -> ag, MaxIterations -> maxiter,
			PrecisionGoal -> pg, WorkingPrecision -> wp}
		];
	FMresult = Apply[FindMinimum, args];

        If[!FreeQ[FMresult, FindMinimum],
                Message[SpatialMedian::findmin];
                Return[$Failed]];
	  
	result = coordvar /. (FMresult[[2]]);

	result
	
   ]

L1[m_, vector_] := Apply[Plus, Map[Sqrt[Apply[Plus, (vector-#)^2]]&, m]]


(* ========= SimplexMedian multivariate robust measure of location ========= *)
Options[SimplexMedian] =
	{AccuracyGoal -> Automatic, 
	 MaxIterations -> 60, PrecisionGoal -> Automatic,
	 WorkingPrecision -> Precision[1.]}


SimplexMedian::invpar =
"Warning: `` `` must be a positive integer.  Default value assumed."

SimplexMedian::findmin =
"Minimization of objective function using FindMinimum failed."

SimplexMedian[m_?MatrixQ, opts___] :=
   Module[{median = iSimplexMedian[m, opts]},
	median /; median =!= $Failed
   ]

iSimplexMedian[m_, opts___] :=
   Module[{n, p, coordvar, coord, unitVec,
	   ag, maxiter, pg, wp, allpsubsets,
	   marginalMedian, volmarg, mean, volmean, init, volinit,
	   savesubvalues, result, detlist, args, 
	   FMresult, vol},
	{n, p} = Dimensions[m];
	coordvar = Table[Unique[coord], {p}]; 
	unitVec = Table[1, {p + 1}];
	{ag, maxiter, pg, wp} =
	    {AccuracyGoal, MaxIterations, PrecisionGoal, WorkingPrecision} /.
		{opts} /. Options[SimplexMedian];
	If[!IntegerQ[maxiter] || maxiter < 1,
	   Message[SimplexMedian::invpar, MaxIterations, maxiter];
	   maxiter = MaxIterations /. Options[SimplexMedian]
        ];
	If[!IntegerQ[wp] || wp < 1,
	   Message[SimplexMedian::invpar, WorkingPrecision, wp];
	   wp = WorkingPrecision /. Options[SimplexMedian]
        ];
	(* For n = 20, p = 4, KSubsets[Range[n], p] has length 4845. *)
	allpsubsets = KSubsets[Range[n], p];

	marginalMedian = Median[m];
	volmarg = totalvolume[m, allpsubsets, marginalMedian];
	mean = Mean[m];
	volmean = totalvolume[m, allpsubsets, mean];
	{init, volinit} = If[volmarg < volmean,
				{marginalMedian, volmarg},
				{mean, volmean}];

	savesubvalues = SubValues[Derivative]; (* save SubValues *)
	Derivative[1][Abs][x_] := Sign[x];

	detlist = Map[Det[ Prepend[Transpose[Append[#, coordvar]],
				 unitVec] ]&,
			Map[m[[#]]&, allpsubsets]  (* p p-vectors *)
		    ];
	detlist = Map[If[Precision[#] < wp,
			   SetPrecision[#, wp],
			   #]&,
		        detlist];
	args = Join[ { Apply[Plus, Map[Abs, detlist]] },
		       Transpose[{coordvar, init}],
		     { AccuracyGoal -> ag, MaxIterations -> maxiter,
			PrecisionGoal -> pg, WorkingPrecision -> wp}
		];
	FMresult = Apply[FindMinimum, args];

        If[!FreeQ[FMresult, FindMinimum],
                Message[SimplexMedian::findmin];
                Return[$Failed]];
	  
	result = coordvar /. (FMresult[[2]]);

	SubValues[Derivative] = savesubvalues; (* restore SubValues *)

	result
	
   ]




(* signed volume of p-dimensional simplex is actually
	signedvolume[m]/Factorial[p] *)

signedvolume[m_] := Det[Prepend[
	Transpose[m], Table[1, {Length[m]}] ]] (* matrix is ((p+1) x p) *)

(* sum of (unsigned) volumes of p-dimensional simplices *)
(* Take the simplex formed by the p-vector "vector" and each set
	of p p-vectors specified in psubsets.  Find the volume of
	that simplex.  Sum the volumes of all simplices so formed. *)
totalvolume[m_, psubsets_, vector_] :=
 With[{unitvec = Table[1, {Length[vector] + 1}]},
  Apply[Plus, Map[
	Abs[Det[
		Prepend[Transpose[ Append[#, vector] ], unitvec ]
	]]&,
        Map[m[[#]]&, psubsets]
  ] ]
 ]

(* KSubsets function from Combinatorica.m *)

KSubsets[l_List,0] := { {} }
KSubsets[l_List,1] := Partition[l,1]
KSubsets[l_List,k_Integer?Positive] := {l} /; (k == Length[l])
KSubsets[l_List,k_Integer?Positive] := {}  /; (k > Length[l])

KSubsets[l_List,k_Integer?Positive] := 
 Join[ Map[(Prepend[#,First[l]])&, KSubsets[Rest[l],k-1]], KSubsets[Rest[l],k] ]



(* ======== ConvexHullMedian multivariate robust measure of location ======= *)
Options[ConvexHullMedian] = {EstimateDOF -> False}

ConvexHullMedian::notimp =
"ConvexHullMedian[matrix] not implemented for matrices n x p, where p > 2."

ConvexHullMedian[coord_?MatrixQ, opts___] :=
   Module[{layer = Layer[coord], innermost, n, median},
     (
      innermost = Last[layer];
      n = Length[innermost];
      median = Apply[Plus, coord[[innermost]]] / n;
      If[TrueQ[EstimateDOF /.
		 {opts} /. Options[ConvexHullMedian]],
	{median, EstimateDOF -> n},
	median
      ]
     ) /; FreeQ[layer, Layer]
   ] /; Length[coord[[1]]] == 2

ConvexHullMedian[m_?MatrixQ, opts___] :=
   Module[{},
	Null /; (Message[ConvexHullMedian::notimp];
		 False)
   ] /; Length[m[[1]]] > 2


(* ================================ Layer ============================== *)
(* 
Layer[{{x11, ... , x1p}, ... , {xn1, ... , xnp}}] gives the indices of
those points lying on the convex layers of the p-dimensional n-point input
set ordered from outermost to innermost.  Layer[{{x11, ... , x1p}, ... ,
{xn1, ... , xnp}}, m] gives the outermost m layers.
*)

Layer[coord_] :=
	Module[{result = iLayer[coord, Infinity]},
	       result /; result =!= $Failed
	]

Layer[coord_, m_Integer?Positive] :=
	Module[{result = iLayer[coord, m]},
			result /; result =!= $Failed
	]

iLayer[coord_, m_] :=
	Module[{n, p, output = {}, i, singlelayer, vertices},
	   {n, p} = Dimensions[coord];
	   vertices = Range[n];
	   If[n <= p+1, Return[{vertices}]];
	   For[i = 1, (i <= m) && (Length[vertices] > p+1), i++,
	      singlelayer = ConvexHull[coord[[vertices]]];
	      If[!FreeQ[singlelayer, ConvexHull],
	         (* do not issue message when Layer is used internally... 
		 Message[Layer::conv, i]; *)
		 Return[$Failed]
	      ];
	      singlelayer = vertices[[singlelayer]];
	      output = Join[output, {singlelayer}];
	      vertices = Complement[vertices, singlelayer];
	   ];
	   If[Length[vertices] > 0 && m == Infinity,
		(* Anywhere from 1 to p+1 points in the interior. *)
		output = Join[output, {vertices}]];
	   output
        ]

Layer::conv = "ConvexHull failed at layer ``."


(* ==== MultivariateTrimmedMean multivariate robust measure of location ===== *)
Options[MultivariateTrimmedMean] = {EstimateDOF -> False}

MultivariateTrimmedMean::notimp =
"MultivariateTrimmedMean[matrix, f] using the convex hull method of trimming
is not implemented for matrices n x p, where p > 2."

MultivariateTrimmedMean::conv = "Unable to find convex layer ``."

MultivariateTrimmedMean::fit = "Unable to interpolate between convex layers."

MultivariateTrimmedMean[m_?MatrixQ, 0, opts___] :=
    If[TrueQ[EstimateDOF /. {opts} /. Options[MultivariateTrimmedMean]],
       {Mean[m], EstimateDOF -> Length[m]},
       Mean[m]
    ]

(* NOTE: need to add support for other MedianMethods, such as SpatialMedian,
	and need to add ability to do convex hulls in 3 or more dimensions.
   Unfortunately, it is tricky to define a multivariate trimmed mean that
   transitions reasonable smoothly from the mean to a multivariate median. *)
MultivariateTrimmedMean[m_?MatrixQ, f_, opts___] :=
 Module[{mean},
   (
   mean
   ) /; 	 (If[Length[m[[1]]] > 2,
		     Message[MultivariateTrimmedMean::notimp];
		     mean = $Failed,
		     (* 2-dimensional data *)
		     mean = iTrimmedHullMean[m, f, opts] 
		   ];
		  mean =!= $Failed)
	
   ] /;  f <= 1 

iTrimmedHullMean[m_, f_, opts___] :=
  Module[{n = Length[m], trimmed = m, level = 1, goalNum, outerNum = 0,
	  layer, currentNum, meano, meani, innerNum, x, fit, meanf},
    goalNum = n f; 
    layer = Layer[trimmed, 1];
    If[ !FreeQ[layer, Layer],
	Message[TrimmedHullMean::conv, level];
	Return[$Failed]];
    layer = layer[[1]];
    currentNum = Length[layer];
    While[goalNum > outerNum + currentNum, 
      (* add most recent layer to outer layers *)
      outerNum += currentNum;
      trimmed = Delete[trimmed, Map[List, layer]];
      (* compute new layer *)
      layer = Layer[trimmed, 1];
      level++;
      If[ !FreeQ[layer, Layer],
	  Message[TrimmedHullMean::conv, level];
	  Return[$Failed]];
      layer = layer[[1]];
      currentNum = Length[layer]
    ];	(* end While *)
    (* interpolate between layers *)
    meano = Apply[Plus, trimmed] / Length[trimmed];
    trimmed = Delete[trimmed, Map[List, layer]];
    If[trimmed === {},
       meanf = meano,
       meani = Apply[Plus, trimmed] / Length[trimmed];
       innerNum = outerNum+currentNum;
       If[innerNum == goalNum,
	  outerNum = innerNum;
	  meanf = meani,
          If[!FreeQ[fit =
	      Fit[{{outerNum, meano}, {innerNum, meani}}, {1, x}, x], Fit],
             Message[TrimmedHullMean::fit];
             Return[$Failed]];
       meanf = fit /. x -> goalNum
       ]
    ];
    If[TrueQ[EstimateDOF /. {opts} /.
		 Options[MultivariateTrimmedMean]],
       (* outerNum points are not involved in the estimate *)
       {meanf, EstimateDOF -> (n-outerNum)},
       (* return only the estimate *)
       meanf
    ]
  ]



(* =========================== PolytopeQuantile ========================== *)
(* ================= multivariate robust measure of location =============== *)
PolytopeQuantile::notimp =
"PolytopeQuantile[matrix, q] not implemented for matrices n x p,
where p > 2."

PolytopeQuantile[m_?MatrixQ, q_] :=
  Module[{layer, innermost, median, n, qlevel, quantile},
   (
   quantile
   ) /; FreeQ[layer  = Layer[m], Layer] &&
	(innermost = Last[layer];
	 median = Apply[Plus, m[[innermost]]] / Length[innermost];
	 If[q == 0,
	    quantile = Polytope[{median}, {1}]; True,
            (* qlevel = {1., (n-n1)/n, (n-n1-n2)/n, ...} corresponding to the 
    	       fact that the 1st layer contains 100% of the points, the 2nd 
	       contains (n-n1)/n of the points, etc. *)
            n = Length[m];
	    qlevel = Drop[N[(n - FoldList[Plus, 0, Map[Length, layer]]) / n],
		 -1];
	    quantile = iPolytopeQuantile[m, layer, innermost, median,
		qlevel, q];
	    quantile =!= $Failed
         ])
  ] /; Length[m[[1]]] == 2 && Head[q]=!=List && NumberQ[N[q]] &&
	FreeQ[N[q], Complex] &&  0<=N[q]<=1 

iPolytopeQuantile[m_, layer_, innermost_, median_, qlevel_, q_] := 
  Module[{n = Length[m], delta, outer, inner, qi, r, quantileVertices,
	  convexhull},
    If[Last[qlevel] - q >= 0,
       If[Length[innermost] == 1,
	  Return[ Polytope[{median}, {1}] ] ]; 
       If[Length[innermost] == 2,
	  delta = (q/Last[qlevel]) Apply[Subtract, m[[ innermost ]]];
	  Return[ Polytope[{median+delta/2, median-delta/2}, {1, 2}] ] ]
    ];
    (* locate quantile between layers for a local fit *)
    outer = 1;
    Scan[If[q > #, Return[outer], outer++]&, Drop[qlevel, 1]];
    (* determine points lying on quantile polygon *)
    Which[outer == Length[layer] || Length[layer[[outer+1]]] == 1,
	     (* quantile corresponds to innermost layer OR
		inner layer composed of one point *)
             qi = If[outer == Length[layer], 0., qlevel[[outer+1]]];
             r = (q-qi)/(qlevel[[outer]]-qi);
             quantileVertices =
	       Map[(median + r (#-median))&, m[[ layer[[outer]] ]] ],
	  True,
	     (* inner layer composed of two or more points *) 
	     Module[{qo, qi, sorted, ls, interquantileSegments,
		     INNER, OUTER, BOTH},
               qo = qlevel[[outer]];
               inner = outer + 1; qi = qlevel[[inner]];
               sorted = PolarSort[Join[
			Map[{INNER, m[[#]]}&, layer[[inner]]],
			Map[{OUTER, m[[#]]}&, layer[[outer]]]
		      ], median];
	       (* Relabel a pair of vertices sharing a common polar angle
		  (one INNER, one OUTER) as BOTH. *)
	       sorted  = sorted //.
		{{x___, {INNER, p1_, a_}, {OUTER, p2_, a_}, y___} :>
		 {x, {BOTH, {p1, p2}, a}, y}};
	       ls = Length[sorted];
               (* The interquantileSegments run between the two convex layers
		  that contain the desired quantile.  Each segment contains a
		  vertex from the INNER layer, a vertex from the OUTER layer, 
		  or vertices from BOTH layers.  The vertices of the desired
		  quantile layer lie on these segments.  *)
	       interquantileSegments = Map[Switch[sorted[[#, 1]],
		        BOTH,
			sorted[[#, 2]],
			INNER,
		        Module[{i=#, j=#, pi, pj, x, y, rule},
		         While[sorted[[i, 1]] === INNER, i = Mod[i, ls]+1 ];
		         While[sorted[[j, 1]] === INNER, j = Mod[j-2, ls]+1 ];
		         pi = If[sorted[[i, 1]]===BOTH,
			   sorted[[i, 2, 2]], sorted[[i, 2]] ];
		         pj = If[sorted[[j, 1]]===BOTH,
			   sorted[[j, 2, 2]], sorted[[j, 2]] ];
	                 rule = Solve[{Det[{{x, y, 1}, Append[pi, 1],
				Append[pj, 1]}] == 0,
		           	     Det[{{x, y, 1}, Append[median, 1],
				Append[sorted[[#, 2]], 1]}] == 0}, {x, y}][[1]];
			 {sorted[[#, 2]], {x, y} /. rule}] (* end Module *),
			OUTER,
		        Module[{i=#, j=#, pi, pj, x, y, rule},
		         While[sorted[[i, 1]] === OUTER, i = Mod[i, ls]+1 ];
		         While[sorted[[j, 1]] === OUTER, j = Mod[j-2, ls]+1 ];
		         pi = If[sorted[[i, 1]]===BOTH,
			   sorted[[i, 2, 1]], sorted[[i, 2]] ];
		         pj = If[sorted[[j, 1]]===BOTH,
			   sorted[[j, 2, 1]], sorted[[j, 2]] ];
	                 rule = Solve[{Det[{{x, y, 1}, Append[pi, 1],
				Append[pj, 1]}] == 0,
		           	     Det[{{x, y, 1}, Append[median, 1],
				Append[sorted[[#, 2]], 1]}] == 0}, {x, y}][[1]];
			 {{x, y} /. rule, sorted[[#, 2]]}] (* end Module *)
                       ]&, Range[ls] ]; (* end Map *)
               (* Now determine the points on the interquantileSegments
		  corresponding to the qth quantile. *)
	       r = (q-qlevel[[inner]])/(qlevel[[outer]]-qlevel[[inner]]);
	       quantileVertices =
	         Map[(#[[1]] + r (#[[2]]-#[[1]]))&, interquantileSegments]
             ] (* end Module *)
    ];	 (* end Which *)
    (* Make quantile locus convex. *)
    If[!FreeQ[convexhull = ConvexHull[quantileVertices], ConvexHull],
	Return[$Failed]
    ];
    Polytope[quantileVertices[[convexhull]], Range[Length[convexhull]]]
  ]


PolytopeQuantile[m_?MatrixQ, q_, opts___] :=
   Module[{},
	   Null /; (Message[PolytopeQuantile::notimp];
		    False)
   ] /; Length[m[[1]]] > 2 && NumberQ[N[q]] && FreeQ[N[q], Complex] &&
	 0<=N[q]<=1

(* Case where q is a list of quantiles. *)
PolytopeQuantile[m_?MatrixQ, q_?VectorQ] :=
  Module[{layer, innermost, median, n, qlevel, quantiles, quantile, scan},
   (
   quantiles
   ) /; FreeQ[layer  = Layer[m], Layer] &&
	(innermost = Last[layer];
	 median = Apply[Plus, m[[innermost]]] / Length[innermost];
         (* qlevel = {1., (n-n1)/n, (n-n1-n2)/n, ...} corresponding to the fact 
    	    that the 1st layer contains 100% of the points, the 2nd contains
	    (n-n1)/n of the points, etc. *)
         n = Length[m];
	 qlevel = Drop[N[(n - FoldList[Plus, 0, Map[Length, layer]]) / n], -1];
	 quantiles = {};
	 scan = Scan[(If[# == 0,
			 quantile = Polytope[{median}, {1}],
			 quantile = iPolytopeQuantile[m, layer, innermost,
				       median, qlevel, #]];
		      If[quantile =!= $Failed,
			 quantiles = Append[quantiles, quantile],
			 Return[$Failed]])&,
		     q];
         scan =!= $Failed)
  ] /; Length[m[[1]]] == 2 && Apply[And, 
	 Map[(NumberQ[#] && FreeQ[#, Complex] && 0<=#<=1)&, N[q]]]


(* =========================== PolytopeQuartiles ========================= *)
(* ================= multivariate robust measure of location =============== *)

PolytopeQuartiles::notimp =
"PolytopeQuartiles[matrix, q] not implemented for matrices n x p,
where p > 2."

PolytopeQuartiles[m_?MatrixQ] :=
	PolytopeQuantile[m, {.25, .50, .75}] /; Length[m[[1]]] == 2

PolytopeQuartiles[m_?MatrixQ] :=
   Module[{},
	   Null /; (Message[PolytopeQuartiles::notimp];
		    False)
   ] /; Length[m[[1]]] > 2


(* ================= computational geometry utilities ==================== *)

PolarSort[l_, median_]:=
   Module[{n=Length[l],p1,p2,in,sorted},
	(* The centroid of the points is interior to the convex hull. *)
	origin=Apply[Plus,Map[#[[2]]&,l]]/n;
	(* 1st component of elements of 'in' is label,
	   2nd component of elements of 'in' is original coordinate,
	   3rd component of elements of 'in' is centered coordinate,
	   4th component of elements of 'in' is polar angle *)
	in = Map[Join[#, {#[[2]]-median, PolarAngle[#[[2]]-median]}]&, l];
	sorted = Sort[in,
                      Function[{p1,p2},
			       p1[[4]]<p2[[4]] ||
			       (* Changed the test p1[[4]]==p2[[4]] to
			       p1[[4]]-p2[[4]]==0 for numerical precision. *)
			      (p1[[4]]-p2[[4]]==0 &&
			      (p1[[3,1]]^2 + p1[[3,2]]^2 <
			       p2[[3,1]]^2 + p2[[3,2]]^2))
		      ] (* end Function *)
	]; (* end Sort *)
	Map[Delete[#, 3]&, sorted]
   ]

PolarAngle[{x_,y_}] := ArcTan[x, y] 

(* left turn if positive, right turn if negative *)
SignOfArea[{x1_,y1_},{x2_,y2_},{x3_,y3_}]:=
  Module[{area = x1(y2-y3) - x2(y1-y3) + x3(y1-y2), prec},
	  If[(prec = Precision[area]) == Infinity || prec == 0,
		  prec = $MachinePrecision];
	  Sign[Chop[area, 10^(1-prec)]]
  ]


(* =========== matrix-valued multivariate dispersion statistics ============ *)

makerectangular[lowertri_] :=
   Module[{p = Length[lowertri], i, j},
	Table[lowertri[[Max[{i, j}], Min[{i, j}]]], {i, p}, {j, p}]
   ]

CovarianceMatrixMLE[m_?MatrixQ] :=
  Module[{mt, mean = Mean[m], p, n, lowertri},
   (
   mt = Transpose[Map[(#-mean)&, m]];
   lowertri = Map[Table[ Apply[Plus, mt[[i]] mt[[#]] ]/n, {i, #}]&, Range[p]];
   makerectangular[lowertri]
   ) /; ({n, p} = Dimensions[m];
	 n > 0)
  ]

Options[CovarianceMatrix] = {ScaleMethod -> StandardDeviation}

CovarianceMatrix::scale = 
"ScaleMethod `` gives zero scale for one or more of the variables."
CovarianceMatrix::cor = 
"ScaleMethod `` unable to give correlation between one or more pairs of
variables."
CovarianceMatrix::notpdef =
"ScaleMethod `` does not yield a positive definite covariance matrix."
CovarianceMatrix::notimp =
"ScaleMethod `` not implemented for covariance matrix between vectors
represented by two data matrices.  Using ScaleMethod -> StandardDeviation."

CovarianceMatrix[m_?MatrixQ, opt___Rule] :=
  Module[{f, cov},
   (
	cov
   ) /; (f = ScaleMethod /. {opt} /. Options[CovarianceMatrix];
	 If[MemberQ[{StandardDeviation, MeanDeviation, MedianDeviation,
                    QuartileDeviation}, f],
	    cov = robustcovmat[m, f],
	    cov = robustcovmat[m, StandardDeviation] ];
	 cov =!= $Failed)
  ]


robustcovmat[m_, StandardDeviation] :=
  Module[{n = Length[m], cov},
	If[n > 1 && 
	   FreeQ[cov = CovarianceMatrixMLE[m], CovarianceMatrixMLE],
		n/(n-1) cov,
		$Failed
	]
  ]

robustcovmat[m_, f_] :=
  Module[{mT = Transpose[m], scale, lowertri, covmat, p = Length[m[[1]]],
		constant = 1/expectedValue[f]^2},
    scale = Map[f, mT];
    Scan[If[!( FreeQ[#, f] && !(# === 0 || # === 0.) ),
	    Message[CovarianceMatrix::scale, f];
	    Return[$Failed]
	 ]&, scale];
    lowertri = Map[Table[
			rc[mT[[i]], mT[[#]], scale[[i]], scale[[#]], f] *
				scale[[i]] scale[[#]] constant ,
			{i, #}]&, Range[p]];
    If[!FreeQ[lowertri, $Failed],
	Message[CovarianceMatrix::cor, f];
	Return[$Failed]
    ];
    covmat = makerectangular[lowertri];
    If[Det[covmat] > 0,
	covmat,
	Message[CovarianceMatrix::notpdef, f];
	$Failed
    ]	
  ]

rc[xlist_, ylist_, scalex_, scaley_, f_] :=
  Module[{scaleplus, scaleminus},
	If[FreeQ[scaleplus = f[xlist/scalex + ylist/scaley], f] &&
	   FreeQ[scaleminus = f[xlist/scalex - ylist/scaley], f],
	    (scaleplus^2 - scaleminus^2)/(scaleplus^2 + scaleminus^2),
	    $Failed
	]
  ]

CovarianceMatrixMLE[m1_?MatrixQ, m2_?MatrixQ] :=
   Module[{mt1, mean1 = Mean[m1], n1, p1, mt2, mean2 = Mean[m2], n2, p2, n},
     (
     mt1 = Transpose[Map[(#-mean1)&, m1]];
     mt2 = Transpose[Map[(#-mean2)&, m2]];
     Table[Apply[Plus, mt1[[i]] mt2[[j]] ]/n, {i, p1}, {j, p2}]
     ) /; ({n1, p1} = Dimensions[m1];
	   {n2, p2} = Dimensions[m2];
	   n1 == n2 && (n = n1; n > 0))
   ]

CovarianceMatrix[m1_?MatrixQ, m2_?MatrixQ, opts___Rule] :=
   Module[{n, f},
     (
	If[(f = (ScaleMethod /. {opts} /. Options[CovarianceMatrix])) =!=
		StandardDeviation,
	   Message[CovarianceMatrix::notimp, f]];
	n/(n-1) cov
     ) /; ((n = Length[m1]) > 1 &&
	   FreeQ[cov = CovarianceMatrixMLE[m1, m2], CovarianceMatrixMLE])
   ]

CovarianceMatrixOfSampleMean[m_?MatrixQ] :=
   Module[{n},
     (
     CovarianceMatrix[m] / n
     ) /; (n = Length[m]) > 1
   ]


(* =========================== DispersionMatrix ======================== *)
(* ====== matrix-valued multivariate robust measure of dispersion ====== *)

Options[DispersionMatrix] = {EstimateDOF -> False}

DispersionMatrix::notimp =
"DispersionMatrix[matrix] not implemented for matrices n x p, where p > 2."
DispersionMatrix::fail =
"Deleting outermost hull leaves only `` points, insufficient for estimating
matrix."
DispersionMatrix::scale =
"Unable to compute scale factor so that determinant of estimate is unbiased
under multinormality."
DispersionMatrix::layer =
"Outer layer of centered and scaled coordinates does not agree with outer
layer of orginal coordinates."

DispersionMatrix[coord_?MatrixQ, opts___] :=
   Module[{layer, innermost, median, outermost, result},
     (
	result
     ) /; FreeQ[layer  = Layer[coord], Layer] &&
	  (innermost = Last[layer];
	   median = Apply[Plus, coord[[innermost]]] / Length[innermost];
	   outermost = First[layer];
	   If[Length[coord]-Length[outermost] <= 2,
	      Message[DispersionMatrix::fail,
		Length[coord]-Length[outermost]];
              False, True] &&
	   FreeQ[result =
		 iDispersionMatrix[coord, median, outermost, opts],
		 $Failed])
   ] /; Length[coord[[1]]] == 2

iDispersionMatrix[coord_, median_, outermost_, opts___] := 
  Module[{outer = Sort[outermost], centered = Map[(#-median)&, coord],
	  n, p, cT, lowertri, s, sinv, svd,
	  u, w, v, sqrtsinv, scaledcentered, newouter, eV},
    {n, p} = Dimensions[coord];
    cT = Transpose[Delete[centered, Map[List, outer]]];
    lowertri = Map[Table[ Apply[Plus, cT[[i]] cT[[#]] ]/(n-Length[outer]),
	{i, #}]&, Range[p]];
    s = makerectangular[lowertri];
    (* s is an initial estimate of the covariance matrix based on all but
	the outermost layer of data *)
    If[!FreeQ[sinv = Inverse[s], Inverse], Return[$Failed]]; 
    If[!FreeQ[svd = SingularValues[sinv], SingularValues], Return[$Failed]];
    {u, w, v} = svd;
    If[!(Dimensions[u] == {p, p} && u == v), Return[$Failed]];
    sqrtsinv = DiagonalMatrix[Sqrt[w]].u;
    (* scale the centered data according to the variances along the
	diagonal of s *)  
    scaledcentered = Map[(sqrtsinv.#)&, centered];  
    If[!FreeQ[newouter = Layer[scaledcentered, 1], Layer],
       Return[$Failed]];
    newouter = Sort[newouter[[1]]];
    (* If the outer layer of the scaled and centered data is not the same
	as the outer layer of the original data, then error.  This is
	a check to make sure the trimming hasn't yielded a covariance matrix 
	that completely distorts the outer layer. *)
    If[newouter =!= outer,
	 Message[DispersionMatrix::layer];
	 Return[$Failed]];
    eV = expectedValueDispersionMatrixRatio[coord, Length[outer]];
    If[eV === $Failed,
       Message[DispersionMatrix::scale];
       eV = 1
    ];
    s *= (eV)^(1/p);
    If[TrueQ[EstimateDOF /. {opts} /. Options[DispersionMatrix]],
	{s, EstimateDOF -> (n-Length[outer])},
        s	
    ]
  ]

expectedValueDispersionMatrixRatio[coord_, m_] :=
  Module[{sTotal = CovarianceMatrix[coord], sinv, detTotal, mean = Mean[coord],
	  scaleOFcoord, trimPos, coordTrimmed, detTrimmed},
	If[!FreeQ[sinv = Inverse[sTotal], Inverse], Return[$Failed]];
	detTotal = Det[sTotal];
	scaleOFcoord = Transpose[{Range[Length[coord]],
				  Map[(#-mean).sinv.(#-mean) &, coord]}];
	(* sort coordinate scales from smallest to largest *)
	scaleOFcoord = Sort[scaleOFcoord, #1[[2]] < #2[[2]] &];
        (* Compute indices of points that give the largest values for scale. *)
        trimPos = Map[Take[#, 1]&, Take[scaleOFcoord, -m]];
	coordTrimmed = Delete[coord, trimPos];
	detTrimmed = Det[CovarianceMatrix[coordTrimmed]];
	detTotal/detTrimmed
  ]

DispersionMatrix[m_?MatrixQ, opts___] :=
   Module[{},
	Null /; (Message[DispersionMatrix::notimp];
		 False)
   ] /; Length[m[[1]]] > 2


(* =========== scalar-valued bivariate dispersion statistics ============ *)

Options[Covariance] = {ScaleMethod -> StandardDeviation}

Covariance[xlist_?VectorQ, ylist_?VectorQ, opt___] :=
   Module[{cov},
      (
      cov
      ) /; (cov = iCovariance[xlist, ylist, opt]) =!= $Failed
   ] /; Length[xlist] == Length[ylist]

iCovariance[xlist_?VectorQ, ylist_?VectorQ, opt___] :=
   Module[{cov, stdx, stdy},
      Switch[ScaleMethod /. {opt} /. Options[Covariance],
	StandardDeviation, pearsoncov[xlist, ylist],
        MeanDeviation, robustcov[xlist, ylist, MeanDeviation],
	MedianDeviation, robustcov[xlist, ylist, MedianDeviation],
	QuartileDeviation, robustcov[xlist, ylist, QuartileDeviation],
	_, pearsoncov[xlist, ylist]
      ]
   ]

pearsoncov[xlist_, ylist_] :=
  Module[{meanx, meany, n = Length[xlist]},
    If[FreeQ[meanx = Mean[xlist], Mean] && FreeQ[meany = Mean[ylist], Mean],
	Apply[Plus, (xlist-meanx) (ylist-meany)] / (n-1),
	$Failed
    ]
  ]

(* scale is adjusted so that if xlist and ylist are bivariate normal, then
	Covariance[xlist, ylist, ScaleMethod -> method] is approximately
	Covariance[xlist, ylist, ScaleMethod -> StandardDeviation]
*)
robustcov[xlist_, ylist_, f_] :=
  Module[{corr, scalex, scaley},
    If[FreeQ[corr = robustcorr[xlist, ylist, f], f] &&
       FreeQ[scalex = f[xlist], f] && FreeQ[scaley = f[ylist], f],
       corr scalex scaley / expectedValue[f]^2,
       $Failed
    ]
  ]

(* 2 NIntegrate[Abs[x]Exp[-x^2/2]/Sqrt[2 Pi], {x, 0, Infinity}] *)
(* expected value of Abs[x] when x is NormalDistribution[0, 1] *)
expectedValue[MeanDeviation] = 0.7978845608018398
(* Quantile[NormalDistribution[0,1], .75] *)
(* 50 % of the deviations are < 0.6744897501960818 and 50 % are less *)
expectedValue[MedianDeviation] = 0.6744897501960818
(* 1/2 the distance between the 25th and 75th quantiles: 0.6744897501960818 *)
expectedValue[QuartileDeviation] = 0.6744897501960818
(* the following is not actually used... *)
expectedValue[StandardDeviation] = 1.


CovarianceMLE[xlist_?VectorQ, ylist_?VectorQ] :=
   Module[{n = Length[xlist], cov},
     (
	cov (n-1)/n	
     ) /; FreeQ[cov = Covariance[xlist, ylist], Covariance]
   ] /; Length[xlist] == Length[ylist]


(* =========== matrix-valued multivariate association statistics ============ *)

Options[CorrelationMatrix] = {ScaleMethod -> StandardDeviation}

CorrelationMatrix::scale =
"ScaleMethod `` gives zero scale for one or more of the variables."
CorrelationMatrix::cor =
"ScaleMethod `` unable to give correlation between one or more pairs of
variables."
CorrelationMatrix::notpdef =
"ScaleMethod `` does not yield a positive definite covariance matrix."
CorrelationMatrix::notimp =
"ScaleMethod `` not implemented for covariance matrix between vectors
represented by two data matrices.  Using ScaleMethod -> StandardDeviation."


CorrelationMatrix[m_?MatrixQ, opt___Rule] := 
 Module[{f, cor},
  (
  cor
  ) /; (f = ScaleMethod /. {opt} /. Options[CorrelationMatrix];
	If[MemberQ[{StandardDeviation, MeanDeviation, MedianDeviation,
		    QuartileDeviation}, f],
	   cor = robustcormat[m, f],
	   cor = robustcormat[m, StandardDeviation] ];
	cor =!= $Failed)
 ]

robustcormat[m_, StandardDeviation] :=
  Module[{mt, mean = Mean[m], sd, p, n, lowertri},
	If[FreeQ[sd = StandardDeviationMLE[m], StandardDeviationMLE] &&
           Apply[And, Map[(!(# === 0 || # === 0.))&, sd]],
             {n, p} = Dimensions[m];
             mt = Transpose[Map[(#-mean)/sd &, m]];
             lowertri = Map[Table[Apply[Plus, mt[[i]] mt[[#]]]/n,
				{i, #}]&, Range[p]];
             makerectangular[lowertri],
             $Failed
        ]        
  ]

robustcormat[m_, f_] :=
  Module[{mT = Transpose[m], scale, lowertri, cormat, p = Length[m[[1]]]},
    scale = Map[f, mT];
    Scan[If[!( FreeQ[#, f] && !(# === 0 || # === 0.) ),
            Message[CorrelationMatrix::scale, f];
            Return[$Failed]
         ]&, scale];
    lowertri = Map[Table[
                        rc[mT[[i]], mT[[#]], scale[[i]], scale[[#]], f],
                        {i, #}]&, Range[p]];
    If[!FreeQ[lowertri, $Failed],
        Message[CorrelationMatrix::cor, f];
        Return[$Failed]
    ];
    cormat = makerectangular[lowertri];
    If[Det[cormat] > 0,
        cormat,
        Message[CorrelationMatrix::notpdef, f];
        $Failed
    ]
  ]



CorrelationMatrix[m1_?MatrixQ, m2_?MatrixQ, opts___Rule] :=
   Module[{mt1, mean1 = Mean[m1], n1, p1, sd1,
	   mt2, mean2 = Mean[m2], n2, p2, sd2, n},
     (
     If[(f = (ScaleMethod /. {opts} /. Options[CorrelationMatrix])) =!=
		StandardDeviation,
	Message[CorrelationMatrix::notimp, f]];
     mt1 = Transpose[Map[(#-mean1)/sd1 &, m1]];
     mt2 = Transpose[Map[(#-mean2)/sd2 &, m2]];
     Table[Apply[Plus, mt1[[i]] mt2[[j]] ]/n, {i, p1}, {j, p2}]
     ) /; ({n1, p1} = Dimensions[m1];
	   {n2, p2} = Dimensions[m2];
	   n1 == n2 && (n = n1; n > 0)) &&
	   FreeQ[sd1 = StandardDeviationMLE[m1], StandardDeviationMLE] &&
	   FreeQ[sd2 = StandardDeviationMLE[m2], StandardDeviationMLE] 
   ]


(* =========================== AssociationMatrix ======================== *)
(* ====== matrix-valued multivariate robust measure of association ====== *)

Options[AssociationMatrix] = {EstimateDOF -> False}

AssociationMatrix::notimp =
"AssociationMatrix[matrix] not implemented for matrices n x p, where p > 2."

AssociationMatrix[coord_?MatrixQ, opts___] :=
   Module[{dof, disp, p = Length[coord[[1]]], returnNum, s, num, scale, r},
     (
     returnNum = TrueQ[dof];
     If[returnNum,
	s = disp[[1]];
	num = EstimateDOF /. disp[[2]],
        s = disp 
     ];
     scale = Map[Sqrt[s[[#, #]]]&, Range[p]];
     r = Table[s[[i, j]]/(scale[[i]] scale[[j]]), {i, p}, {j, p}];
     If[returnNum,
	{r, EstimateDOF -> num},
	r
     ]
     ) /; (dof = EstimateDOF /. {opts} /. Options[AssociationMatrix];
	   FreeQ[disp = DispersionMatrix[coord, EstimateDOF -> dof],
		 DispersionMatrix])
   ] /; Length[coord[[1]]] == 2

AssociationMatrix[m_?MatrixQ, opts___] :=
   Module[{},
	Null /; (Message[AssociationMatrix::notimp];
		 False)
   ] /; Length[m[[1]]] > 2


(* =========== scalar-valued bivariate association statistics ============ *)

Options[Correlation] = {ScaleMethod -> StandardDeviation}

Correlation[xlist_?VectorQ, ylist_?VectorQ, opt___] :=
   Module[{cor},
      (
      cor
      ) /; (cor = iCorrelation[xlist, ylist, opt]) =!= $Failed
   ] /; Length[xlist] == Length[ylist]

iCorrelation[xlist_?VectorQ, ylist_?VectorQ, opt___] :=
   Module[{cov, stdx, stdy},
      Switch[ScaleMethod /. {opt} /. Options[Correlation],
	StandardDeviation, pearsoncorr[xlist, ylist],
        MeanDeviation, robustcorr[xlist, ylist, MeanDeviation],
	MedianDeviation, robustcorr[xlist, ylist, MedianDeviation],
	QuartileDeviation, robustcorr[xlist, ylist, QuartileDeviation],
	_, pearsoncorr[xlist, ylist]
      ]
   ]

pearsoncorr[xlist_, ylist_] :=
  Module[{cov, stdx, stdy},
   If[FreeQ[cov = Covariance[xlist, ylist], Covariance] &&
      FreeQ[stdx = StandardDeviation[xlist], StandardDeviation] &&
      FreeQ[stdy = StandardDeviation[ylist], StandardDeviation],
	cov/(stdx stdy),
	$Failed
   ]
  ]

robustcorr[xlist_, ylist_, f_] :=
  Module[{scalex, scaley, scaleplus, scaleminus},
    If[FreeQ[scalex = f[xlist], f] && !(scalex === 0 || scalex === 0.) &&
       FreeQ[scaley = f[ylist], f] && !(scaley === 0 || scaley === 0.) &&
       FreeQ[scaleplus = f[xlist/scalex + ylist/scaley], f] &&
       FreeQ[scaleminus = f[xlist/scalex - ylist/scaley], f],
     (scaleplus^2 - scaleminus^2)/(scaleplus^2 + scaleminus^2),
     $Failed
    ]
  ]

(* p. 38, Kendall *)
SpearmanRankCorrelation[xlist_?VectorQ, ylist_?VectorQ] :=
  Module[{n, xrank, yrank, xCorrection, yCorrection, SoS},
    (
     {xrank, xCorrection} = rank[xlist];
     {yrank, yCorrection} = rank[ylist];
     SoS = Apply[Plus, (xrank-yrank)^2];
     ( 1/6 (n^3 - n) - SoS - xCorrection - yCorrection ) /
	Sqrt[ (1/6 (n^3 - n) - 2 xCorrection) (1/6 (n^3 - n) - 2 yCorrection) ]
    ) /; (n = Length[xlist]) == Length[ylist] 
  ] 

rank[zlist_] :=
  Module[{z, range = Range[Length[zlist]], tiedrank, repeatCounts, prevVal,
	  ranksOfRepeatsList, correctionTerm, l, meanRank},
    (* Remember original order and sort according to value. *)
    z = Sort[Transpose[{zlist, range}], #1[[1]] < #2[[1]] &];
    (* Associate sorted values with ranks. *)
    z = Transpose[Append[Transpose[z], range]];
    tiedrank = {};  repeatCounts = {};
    prevVal = z[[1, 1]];  ranksOfRepeatsList = {z[[1, 3]]};
    correctionTerm = 0;
    (* Repeated values are assigned a mean rank. *)
    Scan[If[#[[1]] - prevVal == 0,
           (* REPEAT *)
           AppendTo[ranksOfRepeatsList, #[[3]]],
           (* NEW VALUE *)
           l = Length[ranksOfRepeatsList];
           meanRank = Apply[Plus, ranksOfRepeatsList] / l;
           tiedrank = Join[tiedrank, Table[meanRank, {l}]];
           If[l > 1,
	      (* Keep track of the count of this repeat for the correction. *)
	      AppendTo[repeatCounts, l] ];
           prevVal = #[[1]];
           ranksOfRepeatsList = {#[[3]]}
         ]&, Rest[z]];
    l = Length[ranksOfRepeatsList];
    meanRank = Apply[Plus, ranksOfRepeatsList] / l;
    tiedrank = Join[tiedrank,  Table[meanRank, {l}]];
    If[l > 1,  AppendTo[repeatCounts, l] ];
    correctionTerm = Apply[Plus, Map[(#^3 - #)&, repeatCounts]]/12;
    z = Transpose[ Append[Drop[Transpose[z], -1], tiedrank] ];
    (* Reorder according to original order. *)
    z = Sort[z, #1[[2]] < #2[[2]] &];
    {Last[Transpose[z]], correctionTerm}
  ]


KendallRankCorrelation[xlist_?VectorQ, ylist_?VectorQ] :=
  Module[{n, i, j, signs, nc, nd, nx, ny},
    (
	signs = Map[Sign[{Apply[Subtract, xlist[[#]]],
			  Apply[Subtract, ylist[[#]]]}]&,
                    Flatten[Table[{i, j}, {i, n-1}, {j, i+1, n}], 1] ];
        nc = nd = nx = ny = 0;
	Scan[Which[#[[1]]==#[[2]]==0, Null,
		   #[[1]]==#[[2]], nc++,
		   #[[1]]==0, nx++,
		   #[[2]]==0, ny++,
		   True, nd++]&, signs];
	(nc-nd)/Sqrt[(nc+nd+nx)(nc+nd+ny)]
    ) /; (n = Length[xlist]) == Length[ylist]
  ] 




(* =========== scalar-valued multivariate dispersion statistics ============ *)

Options[GeneralizedVariance] = Options[TotalVariation] =  {MLE -> False}

GeneralizedVariance[m_?MatrixQ, opt___] :=
  Module[{c},
    (
    Det[c]
    ) /; If[TrueQ[!(MLE /. {opt} /. Options[GeneralizedVariance])],
		FreeQ[c = CovarianceMatrix[m], CovarianceMatrix],
		FreeQ[c = CovarianceMatrixMLE[m], CovarianceMatrixMLE]]
  ]
	  
TotalVariation[m_?MatrixQ, opt___] :=
  Module[{c},
    (
    Apply[Plus, Map[c[[#, #]]&, Range[Length[ m[[1]] ]]]] 
    ) /; If[TrueQ[!(MLE /. {opt} /. Options[GeneralizedVariance])],
		FreeQ[c = CovarianceMatrix[m], CovarianceMatrix],
		FreeQ[c = CovarianceMatrixMLE[m], CovarianceMatrixMLE]]
  ]

l2[list_] := Sqrt[Apply[Plus, list^2]]
	  
MultivariateMeanDeviation[m_?MatrixQ] :=
  Module[{mean},
    (
    Apply[Plus, Map[l2[#-mean]&, m]] / Length[m]
    ) /; FreeQ[mean = Mean[m], Mean]
  ]

MultivariateMedianDeviation::conv =
"Unable to find median by peeling convex hulls.  Using coordinate-wise median."
MultivariateMedianDeviation::simp =
"Unable to find median by minimizing total simplex volume.  Using
coordinate-wise median."
MultivariateMedianDeviation::spat =
"Unable to find median by minimizing total Euclidean distance.  Using
coordinate-wise median."

Options[MultivariateMedianDeviation] = {MedianMethod -> Median}

MultivariateMedianDeviation[m_?MatrixQ, opt___] :=
  Module[{median, temp},
    median =
      If[Length[m[[1]]] === 1,
	 (* all medians are identical in the univariate case *)
	 Median[m],
	 (* choice of medians in the multivariate case *)
         Switch[MedianMethod /. {opt} /. Options[MedianMethod],
   	  	Median,
		  Median[m],
		ConvexHullMedian,
		  If[FreeQ[temp = ConvexHullMedian[m], 
			   ConvexHullMedian] && temp =!= $Failed,
		     temp,
	  	     Message[MultivariateMedianDeviation::conv];	
		     Median[m]
		  ],
		SimplexMedian,
		  If[FreeQ[temp = SimplexMedian[m], 
			   SimplexMedian] && temp =!= $Failed,
		     temp,
	  	     Message[MultivariateMedianDeviation::simp];	
		     Median[m]
		  ],
		SpatialMedian,
		  If[FreeQ[temp = SpatialMedian[m], 
			   SpatialMedian] && temp =!= $Failed,
		     temp,
	  	     Message[MultivariateMedianDeviation::spat];	
		     Median[m]
		  ],
		_,
		  Median[m]
	 ]
      ];
     Median[Map[l2[#-median]&, m]] 
  ]


(* ======================= ConvexHullArea ======================= *)
(* ====== scalar-valued multivariate robust measure of dispersion ====== *)

ConvexHullArea::notimp =
"ConvexHullArea[matrix] not implemented for matrices n x p,
where p > 2."

ConvexHullArea[coord_?MatrixQ] :=
   Module[{result},
     (
	result
     ) /; FreeQ[hull = ConvexHull[coord], ConvexHull] &&
	   FreeQ[result = iConvexHullArea[coord, hull], $Failed]
   ] /; Length[coord[[1]]] == 2

iConvexHullArea[coord_, hull_] := 
   Module[{mean = Mean[coord], pairs},
     pairs = Partition[Append[hull, First[hull]], 2, 1];
     Apply[Plus,
	   Map[Abs[signedvolume[Append[coord[[#]], mean]]]&,
	       pairs]]/2
   ]


(* ================== multivariate shape statistics ================= *)

CentralMoment[m_?MatrixQ, q_?VectorQ] :=
   Module[{mean = Mean[m]},
	Mean[  Map[Apply[Times, #]&, Map[((#-mean)^q)&, m]]  ]
   ] /; Length[m[[1]]] === Length[q]

MultivariatePearsonSkewness1[m_?MatrixQ] :=
   Module[{s, sinv, mode, mean = Mean[m]},
     (
     (* need factor of 9 to agree with univariate PearsonSkewness1 *)
     9 (mean-mode).sinv.(mean-mode)
     ) /; FreeQ[s = CovarianceMatrix[m], CovarianceMatrix] &&
	  FreeQ[sinv = Inverse[s], Inverse] &&
	  FreeQ[mode = MultivariateMode[m], MultivariateMode] &&
	  mode =!= {}
   ]

MultivariatePearsonSkewness2::conv =
"Unable to find median by peeling convex hulls, using coordinate-wise median."
MultivariatePearsonSkewness2::simp =
"Unable to find median by minimizing total simplex volume, using
coordinate-wise median."
MultivariatePearsonSkewness2::spat =
"Unable to find median by minimizing total Euclidean distance, using
coordinate-wise median."

Options[MultivariatePearsonSkewness2] = {MedianMethod -> Median}

MultivariatePearsonSkewness2[m_?MatrixQ, opts___] :=
   Module[{s, sinv, median, mean = Mean[m], temp},
     (
     median =
       If[Length[m[[1]]]===1,
	  (* all medians are identical in the univariate case *)
	  Median[m],
	  (* choice of medians in multivariate case *)
          Switch[MedianMethod /. {opts} /. Options[MedianMethod],
   	         Median,
		 Median[m],
		 ConvexHullMedian,
		 If[FreeQ[temp = ConvexHullMedian[m], 
			  ConvexHullMedian] && temp =!= $Failed,
		    temp,
	  	    Message[MultivariatePearsonSkewness2::conv];	
		    Median[m]
		 ],
		 SimplexMedian,
		 If[FreeQ[temp = SimplexMedian[m], 
		          SimplexMedian] && temp =!= $Failed,
		    temp,
	  	    Message[MultivariatePearsonSkewness2::simp];	
		    Median[m]
		 ],
		 SpatialMedian,
		 If[FreeQ[temp = SpatialMedian[m], 
		          SpatialMedian] && temp =!= $Failed,
		    temp,
	  	    Message[MultivariatePearsonSkewness2::spat];	
		    Median[m]
		 ],
		 _,
		 Median[m]
          ]
       ];
     (* need factor of 9 to agree with univariate PearsonSkewness2 *)
     9 (mean-median).sinv.(mean-median)
     ) /; FreeQ[s = CovarianceMatrix[m], CovarianceMatrix] &&
	  FreeQ[sinv = Inverse[s], Inverse]
   ]

MultivariateSkewness[m_?MatrixQ] :=
   Module[{s, sinv, mean = Mean[m], m1, n = Length[m]},
     (
     m1 = Map[(#-mean)&, m];
     Apply[Plus, Flatten[Table[If[i==j, 1, 2] (m1[[i]].sinv.m1[[j]])^3,
	{i, n}, {j, i}]]]/n^2
     ) /; FreeQ[s = CovarianceMatrixMLE[m], CovarianceMatrixMLE] &&
          FreeQ[sinv = Inverse[s], Inverse]
   ]

MultivariateKurtosis[m_?MatrixQ] :=
   Module[{s, sinv, mean = Mean[m], m1, n = Length[m]},
     (
     m1 = Map[(#-mean)&, m];
     Apply[Plus, (Map[(#.sinv.#)&, m1])^2]/n
     ) /; FreeQ[s = CovarianceMatrixMLE[m], CovarianceMatrixMLE] &&
	  FreeQ[sinv = Inverse[s], Inverse]
   ]

MultivariateKurtosisExcess[m_?MatrixQ] :=
   Module[{p, kurtosis},
     (
     p = Length[m[[1]]];
     kurtosis - p(p+2)
     ) /; FreeQ[kurtosis = MultivariateKurtosis[m], MultivariateKurtosis]
   ]


(* ========================= data transformations ======================== *)

Options[PrincipalComponents] = {Method -> CovarianceMatrix,
	WorkingPrecision -> $MachinePrecision}

(* Spectral decomposition of symmetric matrix s is
  Transpose[evec].DiagonalMatrix[eval].evec
  where {eval, evec} = Eigensystem[s] .
*)
(* Eigenvectors[CovarianceMatrix[m]] and Eigenvectors[CovarianceMatrixMLE[m]]
	are equivalent because the different scaling factors 
	(i.e., CovarianceMatrix[m] = (n-1)/n CovarianceMatrixMLE[m],
		where n = Length[m[[1]]] )
	are absorbed into Eigenvalues.
   So there is no point in choosing CovarianceMatrix versus CovarianceMatrixMLE.
   The only interesting choices are CovarianceMatrix and CorrelationMatrix.
*)
PrincipalComponents[m_?MatrixQ, opts___] :=
  Module[{method, wp, nm, s, evec, mean},
   (
    (* Change sign of evec so that first principal component is
	positively related to the first original variable.  This
	makes the transformed data easier to interpret. *)
    If[evec[[1, 1]] < 0,
	evec = -evec];
    If[method === CorrelationMatrix,
       standardizedm = Standardize[nm];
       (* Note that Mean[standardizedm] should be approximately {0, ..., 0}. *)
       (* rotate data in standardizednm to principal components;
          Mardia, Kent, Bibby: Sect. 1.5.3, p. 15 *)
       Map[(evec.(#))&, standardizedm],
       (* else *)
       mean = Mean[nm];
       (* rotate data in nm to principal components;
          Mardia, Kent, Bibby: Sect. 1.5.3, p. 15 *)
       Map[(evec.(#-mean))&, nm]
    ]
   )   /;  ({method, wp} = {Method, WorkingPrecision} /. {opts} /.
		 Options[PrincipalComponents];
	    nm = If[Precision[m] === Infinity,
			N[m, wp],
			m];
	    If[method === CovarianceMatrix || method === CovarianceMatrixMLE ||
	       method === CorrelationMatrix,
		FreeQ[s = method[nm], method],
	        Message[PrincipalComponents::meth, method];
		method = CovarianceMatrix;
	        FreeQ[s = method[nm], method]
	    ]) &&
   	    ( FreeQ[evec = Eigenvectors[s], Eigenvectors] )
  ]


PrincipalComponents::meth =
"Possible settings of Method option are CovarianceMatrix or CorrelationMatrix.
Setting Method -> CovarianceMatrix."



(* ====================================================================== *)
End[]

Protect[Mean, GeometricMean, HarmonicMean, RootMeanSquare,
	TrimmedMean, Quantile, InterpolatedQuantile, Median, Mode, Quartiles,
	SampleRange, Variance, VarianceMLE, VarianceOfSampleMean,
	StandardDeviation, StandardDeviationMLE, StandardErrorOfSampleMean,
 	MeanDeviation, MedianDeviation, InterquartileRange,
	QuartileDeviation, CentralMoment,
	Skewness, PearsonSkewness1, PearsonSkewness2,
	QuartileSkewness, Kurtosis, KurtosisExcess];

EndPackage[]


(* :Examples:

Table 1.2.1 "Multivariate Analysis", Mardia, Kent, & Bibby:

data = 
{{77, 82, 67, 67, 81}, {63, 78, 80, 70, 81}, {75, 73, 71, 66, 81},
{55, 72, 63, 70, 68}, {63, 63, 65, 70, 63}, {53, 61, 72, 64, 73},
{51, 67, 65, 65, 68}, {59, 70, 68, 62, 56}, {62, 60, 58, 62, 70},
{64, 72, 60, 62, 45}, {52, 64, 60, 63, 54}, {55, 67, 59, 62, 44},
{50, 50, 64, 55, 63}, {65, 63, 58, 56, 37}, {31, 55, 60, 57, 73},
{60, 64, 56, 54, 40}, {44, 69, 53, 53, 53}, {42, 69, 61, 55, 45},
{62, 46, 61, 57, 45}, {31, 49, 62, 63, 62}, {44, 61, 52, 62, 46},
{49, 41, 61, 49, 64}, {12, 58, 61, 63, 67}, {49, 53, 49, 62, 47},
{54, 49, 56, 47, 53}, {54, 53, 46, 59, 44}, {44, 56, 55, 61, 36},
{18, 44, 50, 57, 81}, {46, 52, 65, 50, 35}, {32, 45, 49, 57, 64},
{30, 69, 50, 52, 45}, {46, 49, 53, 59, 37}, {40, 27, 54, 61, 61},
{31, 42, 48, 54, 68}, {36, 59, 51, 45, 51}, {56, 40, 56, 54, 35},
{46, 56, 57, 49, 32}, {45, 42, 55, 56, 40}, {42, 60, 54, 49, 33},
{40, 63, 53, 54, 25}, {23, 55, 59, 53, 44}, {48, 48, 49, 51, 37},
{41, 63, 49, 46, 34}, {46, 52, 53, 41, 40},
{46, 61, 46, 38, 41}, {40, 57, 51, 52, 31}, {49, 49, 45, 48, 39},
{22, 58, 53, 56, 41}, {35, 60, 47, 54, 33}, {48, 56, 49, 42, 32},
{31, 57, 50, 54, 34}, {17, 53, 57, 43, 51}, {49, 57, 47, 39, 26},
{59, 50, 47, 15, 46}, {37, 56, 49, 28, 45}, {40, 43, 48, 21, 61},
{35, 35, 41, 51, 50}, {38, 44, 54, 47, 24}, {43, 43, 38, 34, 49},
{39, 46, 46, 32, 43}, {62, 44, 36, 22, 42}, {48, 38, 41, 44, 33},
{34, 42, 50, 47, 29}, {18, 51, 40, 56, 30}, {35, 36, 46, 48, 29},
{59, 53, 37, 22, 19}, {41, 41, 43, 30, 33}, {31, 52, 37, 27, 40},
{17, 51, 52, 35, 31}, {34, 30, 50, 47, 36}, {46, 40, 47, 29, 17},
{10, 46, 36, 47, 39}, {46, 37, 45, 15, 30}, {30, 34, 43, 46, 18},
{13, 51, 50, 25, 31}, {49, 50, 38, 23, 9}, {18, 32, 31, 45, 40},
{8, 42, 48, 26, 40}, {23, 38, 36, 48, 15}, {30, 24, 43, 33, 25},
{3, 9, 51, 47, 40}, {7, 51, 43, 17, 22}, {15, 40, 43, 23, 18},
{15, 38, 39, 28, 17}, {5, 30, 44, 36, 18}, {12, 30, 32, 35, 21},
{5, 26, 15, 20, 20}, {0, 40, 21, 9, 14}};


PrincipalComponents[data, Method->CovarianceMatrix]

   PrincipalComponents[,Method->CovarianceMatrix],
	 s = {{305.768, 127.223, 101.579, 106.273, 117.405}, 
    		{127.223, 172.842, 85.1573, 94.6729, 99.012}, 
    		{101.579, 85.1573, 112.886, 112.113, 121.871}, 
    		{106.273, 94.6729, 112.113, 220.38, 155.536}, 
    		{117.405, 99.012, 121.871, 155.536, 297.755}}
	 evals = Eigenvalues[s]
		 {686.99, 202.111, 103.747, 84.6304, 32.1533}
	 sum = Apply[Plus, evals]
		 1109.63
	 evals/sum
		 {0.619115, 0.182142, 0.0934968, 0.0762689, 0.0289766}
	 	(percentages agree with p. 220, Mardia, Kent, & Bibby)

	evecs = Eigenvectors[s]
	evecs = Sign[evecs[[1, 1]]] evecs
		{{0.505446, 0.368349, 0.345661, 0.451123, 0.53465}, 
    		{0.748748, 0.207403, -0.0759081, -0.300888, -0.547782}, 
    		{-0.299789, 0.41559, 0.145318, 0.596626, -0.600276}, 
    		{-0.296184, 0.782888, 0.00323634, -0.51814, 0.175732}, 
    		{0.0793939, 0.188876, -0.92392, 0.285522, 0.151232}}
		
	mean = {38.9545, 50.5909, 50.6023, 46.6818, 42.3068}
	Expand[evecs.({x1, x2, x3, x4, x5}-mean)]
		{-99.4943 + 0.505446 x1 + 0.368349 x2 + 0.345661 x3 +
			 0.451123 x4 + 0.53465 x5,
		1.40223 + 0.748748 x1 + 0.207403 x2 - 0.0759081 x3 - 
			0.300888 x4 - 0.547782 x5, 
		-19.1562 - 0.299789 x1 + 0.41559 x2 + 0.145318 x3 +
			 0.596626 x4 - 0.600276 x5,
		-11.48 - 0.296184 x1 + 0.782888 x2 + 0.00323634 x3 - 
			0.51814 x4 + 0.175732 x5, 
		14.3774 + 0.0793939 x1 + 0.188876 x2 - 0.92392 x3 +
			 0.285522 x4 + 0.151232 x5}
		(agrees with p.218, Mardia, Kent, & Bibby, to within roundoff
			error)


PrincipalComponents[data, Method->CorrelationMatrix]

   PrincipalComponents[,Method->CorrelationMatrix]
	s = {{1., 0.553405, 0.546751, 0.409392, 0.389099}, 
	    {0.553405, 1., 0.609645, 0.485081, 0.436449}, 
	    {0.546751, 0.609645, 1., 0.710806, 0.664736}, 
	    {0.409392, 0.485081, 0.710806, 1., 0.607174}, 
	    {0.389099, 0.436449, 0.664736, 0.607174, 1.}}
	evals = Eigenvalues[s]
		{3.18098, 0.739572, 0.444965, 0.387892, 0.246591}
	sum = Apply[Plus, evals]
		5.
	evals/sum
		{0.636196, 0.147914, 0.088993, 0.0775784, 0.0493182}
	 	(percentages agree with p. 220, Mardia, Kent, & Bibby)

	evecs = Eigenvectors[s]
	evecs = Sign[evecs[[1, 1]]] evecs
		{{0.399605, 0.431419, 0.503282, 0.456994, 0.438244}, 
    		{0.645458, 0.441505, -0.129068, -0.387906, -0.470454}, 
    		{0.620782, -0.705006, -0.037049, -0.136182, 0.312533}, 
    		{-0.145787, 0.298135, -0.108599, -0.666256, 0.658916}, 
    		{-0.130672, -0.181748, 0.846689, -0.422189, -0.234022}}

	(evecs.{u1, u2, u3, u4, u5})[[1]]
		0.399605 u1 + 0.431419 u2 + 0.503282 u3 + 0.456994 u4 + 
			0.438244 u5
		(agrees with bottom of p.219, Mardia, Kent, & Bibby)
	
	sd = StandardDeviation[N[data]]
		{17.4862, 13.1469, 10.6248, 14.8452, 17.2556}

	((evecs.{u1, u2, u3, u4, u5})[[1]]) /.
		 {u1->x1/17.4862, u2->x2/13.1469, u3->x3/10.6248,
			 u4->x4/14.8452, u5->x5/17.2556}

		0.0228526 x1 + 0.0328153 x2 + 0.0473686 x3 + 0.030784 x4 + 
			0.0253972 x5
		(agrees with top of p.220, Mardia, Kent, & Bibby)
		

	





*)


