(*:Mathematica Version: 3.1 *)

(*:Package Version: 1.0 *)

(*:Name: Statistics`MultiDiscreteDistributions` *)

(*:Context: Statistics`MultiDiscreteDistributions` *)

(*:Title: Multivariate Discrete Statistical Distributions *)

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

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

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

(*:Summary:
This package provides properties and functionals of discrete
multivariate probability distributions.
*)

(*:Keywords: discrete distribution, multivariate distribution *)

(*:Requirements: No special system requirements. *)

(*:Warning:
This package extends the definition of several descriptive statistics
functions and the definition of Random.  If the original usage messages
are reloaded, this change will not be reflected in the usage message,
although the extended functionality will remain.
This package also adds rules to Series.
*)

(*:Discussion:
Mean[dist], where dist is the distribution of a random vector, yields a
vector of means, where each element is the mean of the corresponding
element in the random vector.
The same goes for StandardDeviation, Variance, Skewness, Kurtosis,
and KurtosisExcess.
Multivariate functionals that take into account the dependence between
components of a random vector include: CovarianceMatrix,
CorrelationMatrix, MultivariateSkewness, MultivariateKurtosis, and
MultivariateKurtosisExcess.

In this package, k is used to denote the dimensionality of the
multidimensional distribution, as in the Johnson, Kotz, Balakrishnan
reference.
*)

(*:Reference:
  Norman L. Johnson, Samuel Kotz, N. Balakrishnan, Discrete Multivariate
	Distributions, Wiley, 1997.
  Norman L. Johnson, Samuel Kotz, Discrete Distributions, Wiley, 1969.
  Multivariate Analysis, Mardia, Kent, Bibby, Academic Press, 1979.
	(p. 31: definitions of multivariate skewness and kurtosis)
*)

(*:Limitations:
	Quantile, EllipsoidQuantile, and RegionProbability are not
		implemented for MultinomialDistribution,
		NegativeMultinomialDistribution, or MultiPoissonDistribution.

	MultivariateSkewness, MultivariateKurtosis, and ExpectedValue for 
		NegativeMultinomialDistribution are implemented using NSum,
		which does not appear to work very well for anything except
		single sums (sums over one index).   ExpectedValue for
		MultiPoissonDistribution is also implemented using NSum.

	For dimension greater than two, CharacteristicFunction is not
		implemented for MultiPoissonDistribution.

        The vector generalization of BernoulliDistribution,
		MultiBernoulliDistribution, is not implemented.
		See p. 105, Section 36.7, of the Johnson reference.
        The vector generalization of HypergeometricDistribution,
		MultiHypergeometricDistribution, is not implemented.
		See p. 171, Chapter 39, of the Johnson reference.
        The vector generalization of DiscreteUniformDistribution,
		MultiDiscreteUniformDistribution, is not implemented.
		(Note that DiscreteUniformDistribution is a generalized
		HypergeometricDistribution ("Discrete Distributions",
		Johnson & Kotz), so perhaps MultiDiscreteUniformDistribution
		can be obtained as a generalization of
		MultiHypergeometricDistribution.)
	The vector generalization of GeometricDistribution,
		MultiGeometricDistribution, is not implemented.
		See p. 96-97, Section 36.2, of the Johnson reference.
	The vector generalization of LogSeriesDistribution,
		MultiLogSeriesDistribution, is not implemented.
		See p. 157, Section 38.4, of the Johnson reference.

*)

BeginPackage["Statistics`MultiDiscreteDistributions`",
	     (* needed for MultivariateSkewness, MultivariateKurtosis,
		MultivariateKurtosisExcess *)
             "Statistics`MultiDescriptiveStatistics`",
	     (* needed for Mean, Variance, StandardDeviation, Skewness,
		Kurtosis, KurtosisExcess *)
             "Statistics`DescriptiveStatistics`",
	     (* needed for BinomialDistribution *)
	     "Statistics`DiscreteDistributions`",
	     (* needed for CovarianceMatrix and CorrelationMatrix: *)
             "Statistics`Common`MultivariateCommon`",
	     (* needed for Domain, PDF, CDF, CharacteristicFunction,
		RegionProbability: *)
             "Statistics`Common`DistributionsCommon`"
		]

(* vector generalization of BinomialDistribution *)
MultinomialDistribution::usage =
"MultinomialDistribution[n, p] represents the multinomial distribution
for n trials and probability vector p.  A random k-vector {n1, ..., nk}
generated by MultinomialDistribution[n, p], gives the number nj of
occurences of the jth event for the j = 1, ..., k mutually exclusive events
possible at each trial.  The jth element of p is the probability of the
jth event occuring in a given trial.  The index n must be a positive integer,
and the cell probabilities in vector p must sum to unity."

(* vector generalization of NegativeBinomialDistribution *)
NegativeMultinomialDistribution::usage =
"NegativeMultinomialDistribution[n, p] represents the negative multinomial
distribution for success count n and failure probability vector p.  A random 
k-vector {n1, ..., nk} generated by NegativeMultinomialDistribution[n, p] gives 
the number nj of failures of type j that occur before n successes occur.
The jth element of p is the probability of a type j failure in each trial, and
the probability of success in each trial is 1 - Apply[Plus, p]."
 
(* vector generalization of PoissonDistribution (this generalization
	is not unique but is popular)  *)
MultiPoissonDistribution::usage =
"MultiPoissonDistribution[mu0, mu] represents the multiple Poisson
distribution with mean vector mu0 + mu.  A random k-vector {n1, ..., nk}
generated by MultiPoissonDistribution[mu0, mu] gives the number nj of
occurences of the jth event for the j = 1, ..., k mutually exclusive events
possible at each trial.  The mean of nj is mu0 + muj, where muj is the
jth element of mu."





(* ========================================================================
  Extend usage messages for functions defined in Common/DistributionsCommon.m
  Look for the indicated phrase to determine if this has already been done.
*)

If[StringPosition[PDF::usage, "multivariate"] === {},

PDF::usage = PDF::usage <> " " <>
"The dimensionality of x must match the dimensionality of a multivariate
distribution.";
 
CDF::usage = CDF::usage <> " " <>
"The dimensionality of x must match the dimensionality of a multivariate
distribution.";
 
CharacteristicFunction::usage = CharacteristicFunction::usage <> " " <>
"When the distribution is multivariate, the dimensionality of t must match the
dimensionality of the distribution.";
 
RegionProbability::usage = RegionProbability::usage <> " " <>
"RegionProbability[dist, Ellipsoid[mu, radii, dir]] gives the cumulative
probability within Ellipsoid[mu, radii, dir], when the ellipsoid
corresponds to a constant-probability contour of the specified distribution.";

]

(* ========================================================================
  Extend usage messages for functions defined in DescriptiveStatistics.m.
  Look for the indicated phrase to determine if this has already been done.
*)
 
If[StringPosition[Mean::usage, "distribution is multivariate"] === {},

ExpectedValue::usage = ExpectedValue::usage <> " " <>
"ExpectedValue[f, distribution, {x1, x2, ..., xk}] gives the expected value of
the function f of {x1, x2, ..., xk} with respect to the specified multivariate
statistical distribution for {x1, x2, ..., xk}.";
 
Mean::usage = Mean::usage <> " " <>
"Mean[distribution] gives the mean vector when the specified statistical
distribution is multivariate.";

Variance::usage = Variance::usage <> " " <>
"Variance[distribution] gives the diagonal of the covariance matrix when the
specified statistical distribution is multivariate.";

StandardDeviation::usage = StandardDeviation::usage <> " " <>
"StandardDeviation[distribution] gives the square root of the diagonal of the
covariance matrix when the specified statistical distribution is multivariate.";

Skewness::usage = Skewness::usage <> " " <>
"Skewness[distribution] gives the skewness for each of the marginal 
distributions when the specified statistical distribution is multivariate.";

Kurtosis::usage = Kurtosis::usage <> " " <>
"Kurtosis[distribution] gives the kurtosis for each of the marginal 
distributions when the specified statistical distribution is multivariate.";

KurtosisExcess::usage = KurtosisExcess::usage <> " " <>
"KurtosisExcess[distribution] gives the kurtosis excess for each of the 
marginal distributions when the specified statistical distribution is
multivariate.";

]


(* =========================================================================
  Extend usage messages for functions defined in Common/MultivariateCommon.m
  and MultiDescriptiveStatistics.m.  Look for the indicated phrase to
  determine if this has already been done.
*)

If[StringPosition[Mean::usage,
 "specified multivariate statistical distribution"] === {},

CovarianceMatrix::usage = CovarianceMatrix::usage <> " " <>
"CovarianceMatrix[distribution] gives the covariance matrix of the specified
multivariate statistical distribution.";
 
CorrelationMatrix::usage = CorrelationMatrix::usage <> " " <>
"CorrelationMatrix[distribution] gives the correlation matrix of the specified
multivariate statistical distribution.";
 
MultivariateSkewness::usage = MultivariateSkewness::usage <> " " <>
"MultivariateSkewness[distribution] gives the coefficient of skewness of the
specified multivariate statistical distribution.";
 
MultivariateKurtosis::usage = MultivariateKurtosis::usage <> " " <>
"MultivariateKurtosis[distribution] gives the coefficient of kurtosis of the
specified multivariate statistical distribution.";
 
MultivariateKurtosisExcess::usage = MultivariateKurtosisExcess::usage <> " " <>
"MultivariateKurtosisExcess[distribution] gives the kurtosis excess for the
specified multivariate statistical distribution.";

EllipsoidQuantile::usage = EllipsoidQuantile::usage <> " " <>
"EllipsoidQuantile[distribution, q] gives the ellipsoid centered on
Mean[distribution] containing a fraction q of the specified distribution."

]


(* Unprotect MultiDiscreteDistribution.m symbols. *)
Unprotect[MultinomialDistribution]

(* Unprotect DescriptiveStatistics.m symbols. *)
Unprotect[Mean, Variance, StandardDeviation, Skewness, Kurtosis,
         KurtosisExcess]

(* ======================================================================= *)
Begin["`Private`"]

$maxNumOfPart = 5000; (* Controls warning messages about lower bound on
			 domain size. *)
$maxDomainSize = 5000; (* Controls warning messages about domain size. *) 

(* ======================== Multinomial Distribution ====================== *)
(* p. 32, Chap. 35 "Multinomial Distributions",
	"Multivariate Discrete Distributions", Johnson et al *)

Options[MultinomialDistribution] = {}

(* - - - - - - - - -  ParameterQ[MultinomialDistribution[]] - - - - - - - - - *)
MultinomialDistribution/: ParameterQ[
	MultinomialDistribution[n_, p_, nzpos___]] := And[
    If[FreeQ[N[n], Complex], True,
       Message[MultinomialDistribution::realparm, n]; False],
    If[FreeQ[N[p], Complex], True,
       Message[MultinomialDistribution::realparm, p]; False],
    If[N[n] > 0, If[IntegerQ[n], True,
        Message[MultinomialDistribution::posint, n]; False],
        Message[MultinomialDistribution::posint, n]; False, True],
    If[Map[(0 <= N[#] <= 1)&, p] && Apply[Plus, p] == 1, True,
        Message[MultinomialDistribution::probvect, p]; False, True]
]

MultinomialDistribution::realparm = "Parameter `` is expected to be real."

MultinomialDistribution::posint =
"The parameter `1` describing the number of trials is expected to be a
positive integer."

MultinomialDistribution::probvect =
"The parameter `1` is expected to be a vector of cell probabilities,
each entry between 0 and 1, and all entries summing to unity."


(* - - - - - - - - -  DomainQ[MultinomialDistribution[]] - - - - - - - - - - *)
(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)

MultinomialDistribution/: DomainQ[MultinomialDistribution[n_Integer,
	 p_?VectorQ, nzpos___],
	 list_?((MatrixQ[#, IntegerQ])&)] :=
  Module[{v},
    Scan[(v = #;
	  If[!TrueQ[
		Apply[And, Map[(0 <= # <= n)&, v]] && Apply[Plus, v] == n
	     ],
	     Return[False]
	  ])&, list] =!= False
  ] /; Dimensions[list][[2]] == Length[p]

MultinomialDistribution/: DomainQ[MultinomialDistribution[n_Integer,
	 p_?VectorQ, nzpos___],
	  x_?((VectorQ[#, IntegerQ])&)] :=
	Length[x] == Length[p] &&
  TrueQ[Apply[And, Map[(0 <= # <= n)&, x]]] && TrueQ[Apply[Plus, x] == n]




(* - - - -  finding auxiliery parameters of MultinomialDistribution[]] - - - *)
(* NOTE: Make locations of nonzero cell probabilities in
	 the p cell probability vector an auxiliery parameter, so
	that functionals needing those locations (i.e., PDF) can be faster. *)

MultinomialDistribution[n_, p_?VectorQ] :=
	Module[{k = Length[p], zpos, nzpos},
		zpos = Flatten[Position[p, p1_ /; TrueQ[p1 == 0]]];
		nzpos = Complement[Range[k], zpos];
		MultinomialDistribution[n, p, nzpos]
	]

(* Hide nonzero cell probability position vector using Format. *)
MultinomialDistribution /: Format[MultinomialDistribution[n_, p_, nzpos_]] :=
    (
        MultinomialDistribution[n, Short[p]]
    )


(* - - - - - - - - - -  Domain[MultinomialDistribution[]] - - - - - - - - - - *)

MultinomialDistribution/: Domain[MultinomialDistribution[n_Integer,
	 p_?VectorQ, nzpos___]] :=
	Module[{k = Length[p], (* num of possible events *)
	    numOfPart, (* number of partitions of n into at most k integers *)
		part,	(* partitions of n into at most k integers *)
	        estimatedDomainSize, (* size of domain *)
		domain, pos, zerovec},  
	  (* NOTE: Need a way of generating partitions of
		n exclusively of length k or less; lexicographic order
		is not equivalent to ordering according to partition length.
	  *)
	  If[(numOfPart = myPartitionsP[n, k]) > $maxNumOfPart,
	     (* Issue a warning message so users know what they are
			waiting for as the partitions of n are generated. *)	
	     Message[MultinomialDistribution::domszbnd, n, p, numOfPart ]
	  ];
	  If[k/n < .1 && numOfPart < 1000,
	     (* myPartitions1 is based on NextPartition *)
	     part = myPartitions1[n, k, numOfPart],
	     (* myPartitions2 is based on Partitions *)
	     part = myPartitions2[n, k]
	  ];
	  (* --------------------------------------------- *)
	  (* Note: find final size of domain *)
	  part = Map[Join[#, Table[0, {k-Length[#]}]]&, part];
	  estimatedDomainSize = 0;
	  Scan[(estimatedDomainSize +=
			Apply[Multinomial, Map[Length, Split[#]]])&, part];
	  If[estimatedDomainSize > $maxDomainSize,
	     (* Issue a warning message so users know what they are
			waiting for as the domain is generated. *)	
	     Message[MultinomialDistribution::domsz, n, p, estimatedDomainSize ]
	  ];
	  (* --------------------------------------------- *)
	  domain = Flatten[Map[Permutations, part], 1];
	  If[(pos = Position[p, prob_ /; TrueQ[prob == 0]]) === {},
	     domain,
	     pos = Flatten[pos];
	     zerovec = Table[0, {Length[pos]}];	
	     Select[domain, (#[[pos]] === zerovec)&]
	  ]
	]

MultinomialDistribution::domszbnd =
"Warning: the size of the domain of MultinomialDistribution[``, ``] is
bounded below by ``."

MultinomialDistribution::domsz =
"Warning: the size of the domain of MultinomialDistribution[``, ``] is ``."

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
myPartitionsP::usage =
"myPartitionsP[n, k] gives the number of unrestricted partitions of the
integer n into at most k integers.  Equivalently, it gives the number of
unrestricted partitions of the integer n into integers the largest of which is
k."

myPartitionsP[n_, k_] :=
	Module[{j, z},
		Last[ Series[1/Product[(1-z^j), {j, k}], {z, 0, n}][[3]] ]
	]

(* CHECK:
myPartitionsP[10, 5]
	30
Length[Select[Partitions[10], Length[#] <= 5 &]]
	30
*)
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
myPartitions1::usage =
"myPartitions1[n, k] constructs all partitions of integer n of length less
than or equal to k, in reverse lexicographic order."

myPartitions1[n_, k_, numOfPartitions_] :=
 Module[{part = Table[1, {n}], result = {}},
	While[Length[result] < numOfPartitions,
	  part = NextPartition[part];
	  If[Length[part] <= k, AppendTo[result, part]]
	];
	result
 ]

myPartitions2::usage =
"myPartitions2[n, k] constructs all partitions of integer n of length less
than or equal to k, in reverse lexicographic order."

myPartitions2[n_, k_] :=
	Select[Partitions[n], Length[#] <= k &]

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
	

(* Partitions is from DiscreteMath`Combinatorica` by Steven Skiena *)
Partitions::usage = "Partitions[n] constructs all partitions of integer n in
reverse lexicographic order."
Partitions[n_Integer] := Partitions[n,n]
Partitions[n_Integer,_] := {} /; (n<0)
Partitions[n_Integer,1] := { Table[1,{n}] }
Partitions[_,0] := {}
Partitions[n_Integer,maxpart_Integer] :=
        Join[
                Map[(Prepend[#,maxpart])&, Partitions[n-maxpart,maxpart]],
                Partitions[n,maxpart-1]
        ]

(* NextPartition is from DiscreteMath`Combinatorica` by Steven Skiena *)
NextPartition::usage = "NextPartition[p] gives the integer partition following p
in reverse lexicographic order."
NextPartition[p_List] := Join[Drop[p,-1],{Last[p]-1,1}]  /; (Last[p] > 1)
NextPartition[p_List] := {Apply[Plus,p]}  /; (Max[p] == 1)
NextPartition[p_List] :=
        Module[{index,k,m},
                {index} = First[ Position[p,1] ];
                k = p[[index-1]] - 1;
                m = Apply[Plus,Drop[p,index-1]] + k + 1;
                Join[
                        Take[p,index-2],
                        Table[k,{Quotient[m,k]}],
                        If [Mod[m,k] == 0, {}, {Mod[m,k]}]
                ]
        ]


(* - - - - - - - - - -  PDF[MultinomialDistribution[], ] - - - - - - - - - - *)
(* NOTE: do not check parameters on PDF so that PDF will always
        evaluate and be efficient. *)

(* one or more integers out of range, or integers fail to sum to n *)
MultinomialDistribution/: PDF[MultinomialDistribution[n_Integer, p_?VectorQ,
	nzpos___], x_?(VectorQ[#, IntegerQ]&)] :=
		 0 /; !Apply[And, Map[(0 <= # <= n)&, x]] ||
			!(Apply[Plus, x] == n)

(* general case of PDF *)
MultinomialDistribution/: PDF[MultinomialDistribution[n_, p_?VectorQ,
	nzpos_], x_?VectorQ] :=
    ( 
	If[Apply[Plus, x] == n,
	   Evaluate[
	    Apply[Multinomial, x] Apply[Times, (p[[nzpos]])^(x[[nzpos]])]
	   ],
	   0] 	
    ) /; TrueQ[Length[x] == Length[p]]


(* - - - - - - - - - -  CDF[MultinomialDistribution[], ] - - - - - - - - - - *)


(* one or more integers out of range *)
MultinomialDistribution/: CDF[MultinomialDistribution[n_Integer, p_?VectorQ,
	nzpos___], x_?(VectorQ[#, IntegerQ]&), opts___] :=
	0 /; !Apply[And, Map[(0 <= # <= n)&, x]]


(* general case of CDF *)
MultinomialDistribution/: CDF[MultinomialDistribution[n_, p_?VectorQ,
	nzpos___], x_?VectorQ, opts___] :=
  Module[{k = Length[p], y, yvec, prob, arglist, result},
   (
	yvec = Table[Unique[y], {k}];
        prob = PDF[MultinomialDistribution[n, p], yvec];
        arglist = Transpose[{yvec, Table[0, {k}], x}];
        PrependTo[arglist, prob];
        result = Apply[Sum, arglist];
	If[!FreeQ[result, Sum],
		result /. Thread[Rule[yvec, Table[Unique[], {k}]]],
		result]
   )
  ] /; Apply[And, Map[(!NumberQ[#] || (IntegerQ[#] && 0 <= # <= n))&, x]]
        (* NOTE: you don't want the condition "Apply[Plus, x] <= n" here. *)



(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* The marginal distributions of the MultinomialDistribution[n, {p1, ..., pk}]
	are BinomialDistribution[n, pi]. *)

MultinomialDistribution/: Mean[MultinomialDistribution[n_,
	 p_?VectorQ, nzpos___]] :=
	n p /; ParameterQ[MultinomialDistribution[n, p, nzpos]]


(* StandardDeviation yields a vector in multivar case *)
MultinomialDistribution/: StandardDeviation[MultinomialDistribution[
	n_, p_?VectorQ, nzpos___]] :=
   	  Sqrt[n p (1-p)] /; ParameterQ[MultinomialDistribution[n, p, nzpos]]

(* Variance yields a vector in multivar case *)
MultinomialDistribution/: Variance[MultinomialDistribution[
	n_, p_?VectorQ, nzpos___]] := 
		n p (1-p) /; ParameterQ[MultinomialDistribution[n, p, nzpos]]

(* The covariance matrix of the multinomial distribution is singular
	with rank = (Length[p]-1). *)
MultinomialDistribution/: CovarianceMatrix[MultinomialDistribution[
	n_, p_?VectorQ, nzpos___]] :=
		n (DiagonalMatrix[p]-Map[List, p].{p}) /;
				ParameterQ[MultinomialDistribution[n, p, nzpos]]

(* The correlation matrix of the multinomial distribution is singular
	with rank = (Length[p]-1). *)
MultinomialDistribution/: CorrelationMatrix[MultinomialDistribution[
	n_, p_?VectorQ, nzpos___]] :=
   Module[{sdvector, sigma},
	  (
	  DiagonalMatrix[1/sdvector].sigma.DiagonalMatrix[1/sdvector]
	  ) /; FreeQ[sdvector = StandardDeviation[MultinomialDistribution[
		n, p]], StandardDeviation]  &&
	       FreeQ[sigma = CovarianceMatrix[MultinomialDistribution[
		n, p]], CovarianceMatrix]	
   ] /; ParameterQ[MultinomialDistribution[n, p, nzpos]]

MultinomialDistribution/: Skewness[MultinomialDistribution[
	n_, p_?VectorQ, nzpos___]] := (1 - 2 p)/Sqrt[n p (1-p)] /;
				ParameterQ[MultinomialDistribution[n, p, nzpos]]

(* A value for the multivariate coefficient of skewness close to 0
	indicates elliptical symmetry. *)
(* From Mardia: A "suitable population counterpart of the measure of
	multivariate skewness for a random vector x distributed
	MultinormalDistribution[mu, sigma] is the expected value of	
	({x - mu}.Inverse[sigma].Transpose[{y-mu}])^3",  where x and
	y are independent and y is also MultinormalDistribution[mu, sigma]).
*)
MultinomialDistribution/: MultivariateSkewness[MultinomialDistribution[
	n_Integer, p_?VectorQ, nzpos___]] :=
   Module[{inverse,
	   domain,
	   m = Mean[MultinomialDistribution[n, p]],
	   ldom, mskewness, x, y, i, j },
    (
	(* Note that Domain will warn user of a large domain size. *)
	domain = Domain[MultinomialDistribution[n, p]];
	ldom = Length[domain];
	mskewness = 0;
	Do[({x, y} = domain[[{i, j}]];
	    mskewness += ( If[i==j, 1, 2] *
		           PDF[MultinomialDistribution[n, p], x] *
	                   PDF[MultinomialDistribution[n, p], y] *
	    		   (({x-m}.inverse.Transpose[{y-m}])^3)[[1, 1]] )
	   ),	{i, ldom}, {j, i}];
	mskewness
    ) /; FreeQ[inverse = PseudoInverse[CovarianceMatrix[
		MultinomialDistribution[n, p]]], PseudoInverse]
   ] /; ParameterQ[MultinomialDistribution[n, p, nzpos]]

MultinomialDistribution/: Kurtosis[MultinomialDistribution[
	n_, p_?VectorQ, nzpos___]] :=
		 3 + (1 - 6 p (1-p))/(n p (1-p)) /;
				ParameterQ[MultinomialDistribution[n, p, nzpos]]

(* A value of multivariate kurtosis coefficient close to k*(k+2), where k is
the dimensionality of the data, indicates multinormality. *)
(* From Mardia: A "suitable population counterpart of the measure of
	multivariate kurtosis for a random vector x distributed
	MultinormalDistribution[mu, sigma] is the expected value of	
	({x - mu}.Inverse[sigma].Transpose[{x-mu}])^2
*)
MultinomialDistribution/: MultivariateKurtosis[MultinomialDistribution[
	n_Integer, p_?VectorQ, nzpos___]] := 
   Module[{inverse,
	   domain,
	   m = Mean[MultinomialDistribution[n, p]] },
    (
	(* Note that Domain will warn user of a large domain size. *)
        domain = Domain[MultinomialDistribution[n, p]];
	Apply[Plus, Map[ (PDF[MultinomialDistribution[n, p], #] *
		 (({#-m}.inverse.Transpose[{#-m}])^2)[[1, 1]] )&, domain]]
    ) /; FreeQ[inverse = PseudoInverse[CovarianceMatrix[
		MultinomialDistribution[n, p]]], PseudoInverse]
   ] /; ParameterQ[MultinomialDistribution[n, p, nzpos]]
	
MultinomialDistribution/: KurtosisExcess[MultinomialDistribution[
	n_, p_?VectorQ, nzpos___]] :=
		 (1 - 6 p (1-p))/(n p (1-p)) /; 
				ParameterQ[MultinomialDistribution[n, p, nzpos]]

(* A value of multivariate kurtosis excess close to 0 indicates
	multinormality. *) 
MultinomialDistribution/: MultivariateKurtosisExcess[MultinomialDistribution[
	n_Integer, p_?VectorQ, nzpos___]] :=
  Module[{k = Length[p],
	  k1, (* there are only k1 = k-1 independent variates *)
	  multikurtosis},
    (
	k1 = k-1;
	multikurtosis - k1*(k1+2)
    ) /; FreeQ[multikurtosis =
		 MultivariateKurtosis[MultinomialDistribution[n, p]],
		 MultivariateKurtosis]
  ] /; ParameterQ[MultinomialDistribution[n, p, nzpos]]
 
(* NOTE:
  CharacteristicFunction[dist, t] for a random vector x is defined as the
	expected value of Exp[I t.x], where Length[t] == Length[x].
*)
(* See p. 37 of "Discrete Multivariate Distributions". *)
MultinomialDistribution/: CharacteristicFunction[
      MultinomialDistribution[n_, p_?VectorQ, nzpos___], t_?VectorQ] :=
	(
		(p.Exp[I t])^n	
 	) /; TrueQ[Length[t] == Length[p]] &&
		ParameterQ[MultinomialDistribution[n, p, nzpos]]
 

(*  - - - - - - -  Quantile[MultinomialDistribution[], q] - - - - -  - *)
(* NOTE:
Quantile is not implemented for MultinomialDistribution[n, p].
Actually since there is no unique vector x such that
q = CDF[MultinomialDistribution[n, p], x], it is difficult to even
define Quantile[MultinomialDistribution[n, p], q].  One possibility is
the so-called "equicoordinate one-sided percentage points", i.e., the
quantile is {c, c, ..., c} such that CDF[MultinomialDistribution[n, p],
{c, c, ..., c}] = q. *)

MultinomialDistribution/: Quantile[MultinomialDistribution[n_, p_, nzpos___],
   q_] := NotImplemented /; False

(* - - - - -  EllipsoidQuantile[MultinomialDistribution[], ellipsoid] - - - - *)
MultinomialDistribution/: EllipsoidQuantile[MultinomialDistribution[
	n_, p_, nzpos___], q_] := NotImplemented /; False

(* - - - - -  RegionProbability[MultinomialDistribution[], ellipsoid] - - - - *)
MultinomialDistribution/: RegionProbability[MultinomialDistribution[
	n_, p_, nzpos___], EllipsoidQuantile[mu1_, radii_, dir___?MatrixQ]] :=
			NotImplemented /; False

(* - - - - - - -  ExpectedValue[f, MultinomialDistribution[], x] - - - - - - *)

(* 2 arg syntax *)
MultinomialDistribution/: ExpectedValue[f_Function,
  MultinomialDistribution[n_Integer, p_?VectorQ, nzpos___], opts___?OptionQ] :=
  Module[{nn, (* number of arguments of function f *)
	  k = Length[p], (* number of possible events in a trial; dimension
				of MultinomialDistribution[n, p] *)
	  domain},
   (
    Apply[Plus, Map[
		    (Apply[f, #] *
                        PDF[MultinomialDistribution[n, p], #]
		    )&, domain]]
   ) /; (If[Length[f]==1,
	    (* Function with only a body *)
	    nn = Max[Cases[{f}, Slot[z_]->z, Infinity]],
	    (* Function with a list of formal parameters *)
	    nn = Length[f[[1]]]	];
	 nn <= k) &&
	(
	 (* Note that Domain will warn users of a large domain size. *)
	 FreeQ[domain = Domain[MultinomialDistribution[n, p]], Domain]
	)
  ] /; ParameterQ[MultinomialDistribution[n, p, nzpos]]


(* 3 arg syntax *)
MultinomialDistribution/: ExpectedValue[f_,
	MultinomialDistribution[n_Integer, p_?VectorQ, nzpos___], xvec_?VectorQ,
	opts___?OptionQ] :=			
  Module[{domain},
    (
    Apply[Plus, Map[
                    ( (f /. Thread[xvec :> #]) *
                        PDF[MultinomialDistribution[n, p], #]
                    )&, domain]]
    ) /; (
	  (* Note that Domain will warn users of a large domain size. *)
	  FreeQ[domain = Domain[MultinomialDistribution[n, p]], Domain]
	 )
  ] /; Length[xvec] <= Length[p]
	
	
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE: do not check parameters on Random so that Random will always
        evaluate and be efficient. *)

(* single random vector *)
(* Kemp and Kemp (1987) method, p. 68 of "Discrete Multivariate
	Distributions" *)
MultinomialDistribution/: Random[MultinomialDistribution[n_Integer, p_?VectorQ,
		nzpos___]] :=
  Module[{ntot, ptot, result},
	ntot = 0; ptot = 0; result = {};
	Scan[(bvar = Random[BinomialDistribution[n-ntot, #/(1-ptot)]];
	      ntot += bvar;
	      ptot += #;
	      AppendTo[result, bvar])&, Drop[p, -1]];
	Append[result, n-ntot]
  ]	
	     
(* list of random vectors *)
(* NOTE: this RandomArray implementation provides no speed-up tricks such
   as those provided by RandomArray[MultinormalDistribution[]]. *)
MultinomialDistribution/: RandomArray[
	MultinomialDistribution[n_Integer, p_?VectorQ, nzpos___], dim_] :=
  Module[{dim1 = dim, arglist, random},
	If[IntegerQ[dim1], dim1 = {dim1}];
	arglist = Prepend[Map[List, dim1],
			random[MultinomialDistribution[n, p]] ];
	Apply[Table, arglist] /. {random :> Random}
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]



(* ====================================================================== *)
(* ====================================================================== *)
(* ====================================================================== *)
(* ====================================================================== *)
(* ====================================================================== *)

(* =================== Negative  Multinomial Distribution ================= *)
(* p. 93, Chapter 36 "Negative Multinomial and Other Multinomial-Related
	Distributions", "Multivariate Discrete Distributions", Johnson et al *)

Options[NegativeMultinomialDistribution] = {}

(* - - - - - -  ParameterQ[NegativeMultinormalDistribution[]] - - - - - - *)
NegativeMultinomialDistribution/: ParameterQ[
	NegativeMultinomialDistribution[n_, p_, aux___]] := And[
    If[FreeQ[N[n], Complex], True,
       Message[NegativeMultinomialDistribution::realparm, n]; False],
    If[FreeQ[N[p], Complex], True,
       Message[NegativeMultinomialDistribution::realparm, p]; False],
    If[N[n] <= 0, (* NOTE: it is ok for n to be fractional, p. 94 Johnson *)
        Message[NegativeMultinomialDistribution::pos, n]; False,
	True, True],
    If[Map[(0 <= N[#] < 1)&, p] && Apply[Plus, p] < 1, True,
        Message[NegativeMultinomialDistribution::probvect, p]; False, True]
]


NegativeMultinomialDistribution::realparm =
	 "Parameter `` is expected to be real."

NegativeMultinomialDistribution::pos =
"The parameter `1` describing the number of successes is expected to be a
positive, although it need not be an integer." 

NegativeMultinomialDistribution::probvect =
"The parameter `1` is expected to be a vector of failure probabilities,
each entry between 0 and 1, and all entries summing to a quantity
less than unity."


(* - - - - - -  DomainQ[NegativeMultinomialDistribution[]] - - - - - - - *)
(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)

NegativeMultinomialDistribution/: DomainQ[
	NegativeMultinomialDistribution[n_Integer, p_?VectorQ, aux___],
	 list_?((MatrixQ[#, IntegerQ])&)] :=
   Module[{v},
	Scan[(v = #;
	      If[!TrueQ[Apply[And, Map[(0 <= #)&, v]]],
                 Return[False]]
             )&, list] =!= False
   ] /; Dimensions[list][[2]] == Length[p]

NegativeMultinomialDistribution/: DomainQ[
	NegativeMultinomialDistribution[n_Integer, p_?VectorQ, aux___],
	 x_?((VectorQ[#, IntegerQ])&)] :=
		Length[x] == Length[p] && 
		TrueQ[Apply[And, Map[(0 <= #)&, x]]]




(* -  finding auxiliery parameters of NegativeMultinomialDistribution[]] - *)
(* NOTE: Make locations of nonzero cell probabilities in
	 the p cell probability vector an auxiliery parameter, so
	that functionals needing those locations (i.e., PDF) can be faster. *)
NegativeMultinomialDistribution[n_, p_?VectorQ] :=
	Module[{k = Length[p], zpos, nzpos, q},
		zpos = Flatten[Position[p, p1_ /; TrueQ[p1 == 0]]];
		nzpos = Complement[Range[k], zpos];
		q = 1 - Apply[Plus, p]; (* probability of success *)
		NegativeMultinomialDistribution[n, p, q, nzpos]
	]

(* Hide
   success probability q and nonzero cell probability position vector nzpos
   using Format. *)
NegativeMultinomialDistribution /: Format[
	NegativeMultinomialDistribution[n_, p_, q_, nzpos_]] :=
    (
        NegativeMultinomialDistribution[n, Short[p]]
    )


(*  - - - - - - -  Domain[NegativeMultinomialDistribution[]] - - - - - - - *)

NegativeMultinomialDistribution/: Domain[
	NegativeMultinomialDistribution[n_Integer,
	 p_?VectorQ, aux___]] := 
        Module[{k = Length[p], domain, range = (Head[Range::range] === $Off)},
          Off[Range::range];
          domain = Table[Range[0, Infinity], {k}];
          If[!range, On[Range::range]];
          domain
        ] /; ParameterQ[NegativeMultinomialDistribution[n, p]]


(* - - - - - - -  PDF[NegativeMultinomialDistribution[], ] - - - - - - - *)
(* NOTE: do not check parameters on PDF so that PDF will always
        evaluate and be efficient. *)

(* out of range or nonintegral *)
NegativeMultinomialDistribution/:
	 PDF[NegativeMultinomialDistribution[n_Integer, p_?VectorQ,
	aux___], x_] := 0 /;
        ParameterQ[NegativeMultinomialDistribution[n, p, aux]] &&
	Apply[And, Map[
		 ( NumberQ[N[#]] && (!IntegerQ[#] || Negative[#]) )&, x]]

(* general case of PDF *)
NegativeMultinomialDistribution/:
	 PDF[NegativeMultinomialDistribution[n_, p_?VectorQ,
	q_, nzpos_], x_?VectorQ] :=
	  (
		Apply[Multinomial, Join[x, {n-1}]] q^n *
			Apply[Times, (p[[nzpos]])^(x[[nzpos]])] 	
	  ) /; TrueQ[Length[x] == Length[p]]


(* - - - - - - -  CDF[NegativeMultinomialDistribution[], ] - - - - - - - *)

(* one or more integers out of range *)
NegativeMultinomialDistribution/: CDF[NegativeMultinomialDistribution[
	n_Integer, p_?VectorQ, aux___],
          x_?(VectorQ[#, IntegerQ]&), opts___] := 
	0 /; !Apply[And, Map[(0 <= #)&, x]]


(* general case of CDF *)
NegativeMultinomialDistribution/:
	 CDF[NegativeMultinomialDistribution[n_, p_?VectorQ, aux___],
	  x_?VectorQ, opts___] := 
  Module[{k = Length[p], y, yvec, prob, arglist, result},
   (
	yvec = Table[Unique[y], {k}];
        prob = PDF[NegativeMultinomialDistribution[n, p], yvec];
        arglist = Transpose[{yvec, Table[0, {k}], x}];
        PrependTo[arglist, prob];
	result = Apply[Sum, arglist];
	If[!FreeQ[result, Sum],
                result /. Thread[Rule[yvec, Table[Unique[], {k}]]],
                result]
   )
  ] /; Apply[And, Map[(!NumberQ[#] || (IntegerQ[#] && 0 <= #))&, x]]



(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* The marginal distributions of the
	 NegativeMultinomialDistribution[n, {p1, ..., pk}]
	are NegativeBinomialDistribution[n, pi/q] 
	where q = 1-(p1+...+pk) is the probability of success in a given trial.
	(p. 100, Johnson)
*)

NegativeMultinomialDistribution/: Mean[
	NegativeMultinomialDistribution[n_,
	 p_?VectorQ, q_, aux___]] := n p/q /; 
		ParameterQ[NegativeMultinomialDistribution[n, p, q, aux]]

(* StandardDeviation yields a vector in multivar case *)
NegativeMultinomialDistribution/: StandardDeviation[
	NegativeMultinomialDistribution[n_, p_?VectorQ, q_, aux___]] :=
   		Sqrt[n p (q + p)]/q /;
		ParameterQ[NegativeMultinomialDistribution[n, p, q, aux]]

(* Variance yields a vector in multivar case *)
(* Using Eq. (36.27) for factorial moments and relationship between
	central and factorial moments. *)
NegativeMultinomialDistribution/: Variance[NegativeMultinomialDistribution[
	n_, p_?VectorQ, q_, aux___]] := 
		n p (q + p)/q^2 /;
		ParameterQ[NegativeMultinomialDistribution[n, p, q, aux]]

NegativeMultinomialDistribution/: CovarianceMatrix[
	NegativeMultinomialDistribution[
	n_, p_?VectorQ, q_, aux___]] :=
		n/q^2 (q DiagonalMatrix[p]+Map[List, p].{p}) /;
		ParameterQ[NegativeMultinomialDistribution[n, p, q, aux]]

NegativeMultinomialDistribution/:
	 CorrelationMatrix[NegativeMultinomialDistribution[
	n_, p_?VectorQ, aux___]] :=
   Module[{sdvector, sigma},
	  (
	  DiagonalMatrix[1/sdvector].sigma.DiagonalMatrix[1/sdvector]
	  ) /; FreeQ[sdvector = StandardDeviation[
			NegativeMultinomialDistribution[
		n, p, aux]], StandardDeviation]  &&
	       FreeQ[sigma = CovarianceMatrix[
			NegativeMultinomialDistribution[
		n, p, aux]], CovarianceMatrix]	
   ] /; ParameterQ[NegativeMultinomialDistribution[n, p, aux]]

(* Using Eq. (36.27) for factorial moments, relationship between
        central and factorial moments, and definition of skewness in
	terms of central moments. *)
NegativeMultinomialDistribution/: Skewness[NegativeMultinomialDistribution[
	n_, p_?VectorQ, q_, aux___]] :=
		(2*p + q)/Sqrt[n*p*(p + q)] /;
		ParameterQ[NegativeMultinomialDistribution[n, p, q, aux]]

(* - - - - - - - - - - - - -  MultivariateSkewness - - - - - - - - - - - - - *)
If[FreeQ[Options[MultivariateSkewness], NSumTerms],
	(* add Options[NSum] *)
   Options[MultivariateSkewness] = Join[Options[MultivariateSkewness],
					Options[NSum]]
	]

(* A value for the multivariate coefficient of skewness close to 0
	indicates elliptical symmetry. *)
(* From Mardia: A "suitable population counterpart of the measure of
	multivariate skewness for a random vector x distributed
	MultinormalDistribution[mu, sigma] is the expected value of	
	({x - mu}.Inverse[sigma].Transpose[{y-mu}])^3",  where x and
	y are independent and y is also MultinormalDistribution[mu, sigma]).
*)
NegativeMultinomialDistribution/: MultivariateSkewness[
	NegativeMultinomialDistribution[
	n_Integer, p_?VectorQ, aux___], opts___] :=
   Module[{inverse, k = Length[p], y, yvec, x, xvec, summand,
	   zerovec, infinityvec, arglist,
	   m = Mean[NegativeMultinomialDistribution[n, p, aux]],
	   accgoal, compiled, method, extraterms, terms, precgoal,
           verify, wprec, wynn},
    (
	xvec = Table[Unique[x], {k}];
	yvec = Table[Unique[y], {k}];
        summand = PDF[NegativeMultinomialDistribution[n, p], xvec] *
		  PDF[NegativeMultinomialDistribution[n, p], yvec] *
		  (({xvec-m}.inverse.Transpose[{yvec-m}])^3)[[1, 1]];
	zerovec = Table[0, {k}]; infinityvec = Table[Infinity, {k}];
	arglist = Join[Transpose[{xvec, zerovec, infinityvec}],
			Transpose[{yvec, zerovec, infinityvec}]];
	PrependTo[arglist, summand];
	{accgoal, compiled, method, extraterms, terms, precgoal,
         verify, wprec, wynn} = {AccuracyGoal, Compiled, Method,
                NSumExtraTerms, NSumTerms, PrecisionGoal,
                VerifyConvergence, WorkingPrecision, WynnDegree} /. {opts} /.
                        Options[MultivariateKurtosis];
        Apply[NSum, Join[arglist,
                {AccuracyGoal->accgoal, Compiled->compiled,
                Method -> method, NSumExtraTerms -> extraterms,
                NSumTerms -> terms, PrecisionGoal -> precgoal,
                VerifyConvergence -> verify, WorkingPrecision -> wprec,
                WynnDegree -> wynn}]
	]
    ) /; FreeQ[inverse = Inverse[CovarianceMatrix[
		NegativeMultinomialDistribution[n, p, aux]]], Inverse]
   ] /; ParameterQ[NegativeMultinomialDistribution[n, p, aux]] 
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)


(* Using Eq. (36.27) for factorial moments, relationship between
        central and factorial moments, and definition of skewness in
	terms of central moments. *)
NegativeMultinomialDistribution/: Kurtosis[NegativeMultinomialDistribution[
	n_, p_?VectorQ, q_, aux___]] := 
		3 + (6 p^2 + 6 p q + q^2)/(n p (p+q)) /;
		ParameterQ[NegativeMultinomialDistribution[n, p, q, aux]]

(* - - - - - - - - - - - - - MultivariateKurtosis - - - - - - - - - - - - - - *)
If[FreeQ[Options[MultivariateKurtosis], NSumTerms],
	(* add Options[NSum] *)
   Options[MultivariateKurtosis] = Join[Options[MultivariateKurtosis],
					Options[NSum]]
   ]

(* A value of multivariate kurtosis coefficient close to k*(k+2), where k is
the dimensionality of the data, indicates multinormality. *)
(* From Mardia: A "suitable population counterpart of the measure of
	multivariate kurtosis for a random vector x distributed
	MultinormalDistribution[mu, sigma] is the expected value of	
	({x - mu}.Inverse[sigma].Transpose[{x-mu}])^2
*)
NegativeMultinomialDistribution/: MultivariateKurtosis[
	NegativeMultinomialDistribution[
	n_Integer, p_?VectorQ, aux___], opts___?OptionQ] := 
   Module[{inverse, k = Length[p], y, yvec, summand, arglist,
	   m = Mean[NegativeMultinomialDistribution[n, p, aux]],
	   accgoal, compiled, method, extraterms, terms, precgoal,
	   verify, wprec, wynn },
    (
	yvec = Table[Unique[y], {k}];
	summand = PDF[NegativeMultinomialDistribution[n, p, aux], yvec] *
		(({yvec-m}.inverse.Transpose[{yvec-m}])^2)[[1, 1]];
	arglist = Map[{#, 0, Infinity}&, yvec];
	PrependTo[arglist, summand];
	{accgoal, compiled, method, extraterms, terms, precgoal,
	 verify, wprec, wynn} = {AccuracyGoal, Compiled, Method,
		NSumExtraTerms, NSumTerms, PrecisionGoal,
		VerifyConvergence, WorkingPrecision, WynnDegree} /. {opts} /.
			Options[MultivariateKurtosis];
	Apply[NSum, Join[arglist,
		{AccuracyGoal->accgoal, Compiled->compiled, 
		Method -> method, NSumExtraTerms -> extraterms,
		NSumTerms -> terms, PrecisionGoal -> precgoal,
		VerifyConvergence -> verify, WorkingPrecision -> wprec,
		WynnDegree -> wynn}] 
	]
    ) /; FreeQ[inverse = Inverse[CovarianceMatrix[
		NegativeMultinomialDistribution[n, p, aux]]], Inverse]
   ] /; ParameterQ[NegativeMultinomialDistribution[n, p, aux]] 
	
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)

NegativeMultinomialDistribution/: KurtosisExcess[
	NegativeMultinomialDistribution[
	n_, p_?VectorQ, q_, aux___]] := 
		(6 p^2 + 6 p q + q^2)/(n p (p+q))

(* - - - - - - - - - - - - - MultivariateKurtosisExcess -  - - - - - - - - - *)
If[FreeQ[Options[MultivariateKurtosisExcess], NSumTerms],
	(* add Options[NSum] *)
   Options[MultivariateKurtosisExcess] =
		 Join[Options[MultivariateKurtosisExcess],
					Options[NSum]]
   ]

(* A value of multivariate kurtosis excess close to 0 indicates
	multinormality. *) 
NegativeMultinomialDistribution/: MultivariateKurtosisExcess[
	NegativeMultinomialDistribution[
	n_Integer, p_?VectorQ, aux___], opts___?OptionQ] :=
  Module[{k = Length[p], (* there are k independent variates *)
	  multikurtosis},
    (
	multikurtosis - k*(k+2)
    ) /; FreeQ[multikurtosis =
	 MultivariateKurtosis[NegativeMultinomialDistribution[n, p, aux], opts],
		 MultivariateKurtosis]
  ] 
 
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE:
  CharacteristicFunction[dist, t] for a random vector x is defined as the
	expected value of Exp[I t.x], where Length[t] == Length[x].
*)
(* See p. 101 of "Discrete Multivariate Distributions". *)
NegativeMultinomialDistribution/: CharacteristicFunction[
      NegativeMultinomialDistribution[n_,
	 p_?VectorQ, q_, aux___], t_?VectorQ] :=
	(
		((1 - p.Exp[I t])/q)^(-n)
 	) /; TrueQ[Length[t] == Length[p]]


(*  - - - -  Quantile[NegativeMultinomialDistribution[], q] - - - - - - *)
(* NOTE:
Quantile is not implemented for NegativeMultinomialDistribution[n, p].
Actually since there is no unique vector x such that
q = CDF[NegativeMultinomialDistribution[n, p], x], it is difficult to even
define Quantile[NegativeMultinomialDistribution[n, p], q].  One possibility is
the so-called "equicoordinate one-sided percentage points", i.e., the
quantile is {c, c, ..., c} such that CDF[NegativeMultinomialDistribution[n, p],
{c, c, ..., c}] = q. *)

NegativeMultinomialDistribution/: Quantile[
	NegativeMultinomialDistribution[n_, p_, aux___], q_] :=
   		NotImplemented /; False

(* - -  EllipsoidQuantile[NegativeMultinomialDistribution[], ellipsoid] - - *)
NegativeMultinomialDistribution/: EllipsoidQuantile[
	NegativeMultinomialDistribution[
	n_, p_, aux___], q_] := NotImplemented /; False

(* - -  RegionProbability[NegativeMultinomialDistribution[], ellipsoid] - - *)
NegativeMultinomialDistribution/: RegionProbability[
	NegativeMultinomialDistribution[
	n_, p_, aux___], EllipsoidQuantile[mu1_, radii_, dir___?MatrixQ]] :=
			NotImplemented /; False

(* - - - -  ExpectedValue[f, NegativeMultinomialDistribution[], x] - - - - - *)
If[FreeQ[Options[ExpectedValue], NSumTerms],
	(* add Options[NSum] *)
   Options[ExpectedValue] = Join[Options[ExpectedValue], Options[NSum]]
]

(* 2 arg syntax *)
NegativeMultinomialDistribution/: ExpectedValue[f_Function,
  NegativeMultinomialDistribution[n_Integer,
	 p_?VectorQ, aux___], opts___?OptionQ] :=
  Module[{nn, (* number of arguments of function f *)
	  k = Length[p], (* number of possible events in a trial; dimension
				of NegativeMultinomialDistribution[n, p] *)
	  y, yvec, summand, arglist,
	  accgoal, compiled, method, extraterms, terms, precgoal,
         verify, wprec, wynn
	 },
   (
    yvec = Table[Unique[y], {k}];
    summand = Apply[f, yvec] * PDF[NegativeMultinomialDistribution[
	n, p, aux], yvec];
    arglist = Map[{#, 0, Infinity}&, yvec];
    PrependTo[arglist, summand];
    {accgoal, compiled, method, extraterms, terms, precgoal,
         verify, wprec, wynn} = {AccuracyGoal, Compiled, Method,
                NSumExtraTerms, NSumTerms, PrecisionGoal,
                VerifyConvergence, WorkingPrecision, WynnDegree} /. {opts} /.
                        Options[ExpectedValue];
    Apply[NSum, Join[arglist,
                {AccuracyGoal->accgoal, Compiled->compiled,
                Method -> method, NSumExtraTerms -> extraterms,
                NSumTerms -> terms, PrecisionGoal -> precgoal,
                VerifyConvergence -> verify, WorkingPrecision -> wprec,
                WynnDegree -> wynn}]
	]
   ) /; (If[Length[f]==1,
	    (* Function with only a body *)
	    nn = Max[Cases[{f}, Slot[z_]->z, Infinity]],
	    (* Function with a list of formal parameters *)
	    nn = Length[f[[1]]]	];
	 nn <= k)
  ] /; ParameterQ[NegativeMultinomialDistribution[n, p, aux]]


(* 3 arg syntax *)
NegativeMultinomialDistribution/: ExpectedValue[f_,
	NegativeMultinomialDistribution[n_Integer,
		 p_?VectorQ, aux___], xvec_?VectorQ,
	opts___?OptionQ] :=			
  Module[{summand, arglist,
	  accgoal, compiled, method, extraterms, terms, precgoal,
         verify, wprec, wynn},
    summand = f * PDF[NegativeMultinomialDistribution[
	n, p, aux], xvec];
    arglist = Map[{#, 0, Infinity}&, xvec];
    PrependTo[arglist, summand];
    {accgoal, compiled, method, extraterms, terms, precgoal,
         verify, wprec, wynn} = {AccuracyGoal, Compiled, Method,
                NSumExtraTerms, NSumTerms, PrecisionGoal,
                VerifyConvergence, WorkingPrecision, WynnDegree} /. {opts} /.
                        Options[ExpectedValue];
    Apply[NSum, Join[arglist,
                {AccuracyGoal->accgoal, Compiled->compiled,
                Method -> method, NSumExtraTerms -> extraterms,
                NSumTerms -> terms, PrecisionGoal -> precgoal,
                VerifyConvergence -> verify, WorkingPrecision -> wprec,
                WynnDegree -> wynn}]
	]
  ] /; Length[xvec] <= Length[p] &&
	ParameterQ[NegativeMultinomialDistribution[n, p, aux]]
	
	
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE: do not check parameters on Random so that Random will always
        evaluate and be efficient. *)

negativemultinomial = Compile[{{r, _Real}, {p, _Real, 1}},
        Module[{nsucc = 0, nfail = Table[0, {Length[p]}], p1,
		cumprob = FoldList[Plus, 0, p], psum},
	  psum = Last[cumprob];
	  p1 = Random[];
	  If[p1 > psum,
	     nsucc++,
	     nfail[[ Position[Sort[Prepend[cumprob, p1]], p1][[1, 1]]-1 ]] ++
	  ];
	  While[nsucc < r, 
	    p1 = Random[];
            If[p1 > psum,
               nsucc++,
               nfail[[ Position[Sort[Prepend[cumprob, p1]], p1][[1, 1]]-1 ]] ++
            ]
	  ];
	  nfail
	]	]
	

(* single random vector *)
NegativeMultinomialDistribution/: Random[
	NegativeMultinomialDistribution[n_Integer, p_?VectorQ, aux___]] :=
		negativemultinomial[n, p]
	     
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE: do not check parameters on Random so that Random will always
        evaluate and be efficient. *)

negativemultinomial2 = Compile[{{r, _Real}, {p, _Real, 1},
			{cumprob, _Real, 1}, {psum, _Real}},
        Module[{nsucc = 0, nfail = Table[0, {Length[p]}], p1},
	  p1 = Random[];
	  If[p1 > psum,
	     nsucc++,
	     nfail[[ Position[Sort[Prepend[cumprob, p1]], p1][[1, 1]]-1 ]] ++
	  ];
	  While[nsucc < r, 
	    p1 = Random[];
            If[p1 > psum,
               nsucc++,
               nfail[[ Position[Sort[Prepend[cumprob, p1]], p1][[1, 1]]-1 ]] ++
            ]
	  ];
	  nfail
	]	]

(* list of random vectors *)
NegativeMultinomialDistribution/: RandomArray[
	NegativeMultinomialDistribution[n_Integer,
		 p_?VectorQ, q_, aux___], dim_] :=
  Module[{dim1 = dim, cumprob = FoldList[Plus, 0, p], psum,
		negativemultinomial3},
	psum = Last[cumprob];
	If[IntegerQ[dim1], dim1 = {dim1}];
        arglist = Prepend[Map[List, dim1],
		negativemultinomial3[n, p, cumprob, psum] ];
	Apply[Table, arglist] /. {negativemultinomial3 :> negativemultinomial2}
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]



(* ====================================================================== *)
(* ====================================================================== *)
(* ====================================================================== *)
(* ====================================================================== *)

(* ====================== Multiple Poisson Distribution ================= *)
(* p. 124, Chap. 37 "Multivariate Poisson Distributions",
	"Multivariate Discrete Distributions", Johnson et al *)

Options[MultiPoissonDistribution] = {}

(* - - - - - - - - -  ParameterQ[MultiPoissonDistribution[]] - - - - - - - - *)
MultiPoissonDistribution/: ParameterQ[
	MultiPoissonDistribution[mu0_, mu_, aux___]] := And[
    If[FreeQ[N[mu0], Complex], True,
       Message[MultiPoissonDistribution::realparm, mu0]; False],
    If[FreeQ[N[mu], Complex], True,
       Message[MultiPoissonDistribution::realparm, mu]; False],
    If[N[mu0] <= 0, 
	Message[MultiPoissonDistribution::posparm, mu0]; False,
	True, True],
    If[Apply[And, Map[(N[#] > 0)&, mu]],  True,
        Message[MultiPoissonDistribution::posparm, mu0]; False, True]
]

MultiPoissonDistribution::realparm = "Parameter `` is expected to be real."

MultiPoissonDistribution::posparm =
"Parameter `1` is expected to be positive."



(* - - - - - - - - -  DomainQ[MultiPoissonDistribution[]] - - - - - - - - - - *)
(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)

MultiPoissonDistribution/: DomainQ[
        MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___],
         list_?((MatrixQ[#, IntegerQ])&)] :=
   Module[{v},
	Scan[(v = #;
              If[!TrueQ[Apply[And, Map[(0 <= #)&, v]]],
                 Return[False]]
             )&, list] =!= False
   ] /; Dimensions[list][[2]] == Length[mu]


MultiPoissonDistribution/: DomainQ[
        MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___],
	x_?((VectorQ[#, IntegerQ])&)] :=
                Length[x] == Length[mu] &&
                TrueQ[Apply[And, Map[(0 <= #)&, x]]]


(* - - - -  finding auxiliery parameters of MultiPoissonDistribution[]] - - - *)
(* NOTE: to date haven't found a need for auxiliery parameters for this
	distribution *)

MultiPoissonDistribution /: Format[
	MultiPoissonDistribution[mu0_, mu_?VectorQ, aux__]] :=
    (
        MultiPoissonDistribution[mu0, Short[mu]]
    )



(* - - - - - - - - - -  Domain[MultiPoissonDistribution[]] - - - - - - - - - *)

MultiPoissonDistribution/: Domain[
        MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___]] :=
        Module[{k = Length[mu], domain, range = (Head[Range::range] === $Off)},
          Off[Range::range];
          domain = Table[Range[0, Infinity], {k}];
          If[!range, On[Range::range]];
          domain
        ] /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]



(* - - - - - - - - - -  PDF[MultiPoissonDistribution[], ] - - - - - - - - - - *)
(* NOTE: do not check parameters on PDF so that PDF will always
        evaluate and be efficient. *)

(* one or more negative integers *)
MultiPoissonDistribution/: PDF[MultiPoissonDistribution[mu0_, mu_?VectorQ,
	aux___], x_?((VectorQ[#, IntegerQ])&)] := 0 /;
		Apply[Or, Map[(# < 0)&, x]]

(* univariate Poisson *)
MultiPoissonDistribution/: PDF[MultiPoissonDistribution[mu0_, mu_?VectorQ,
	aux___], x_?VectorQ] :=
	 Module[{mu1, x1},  
		mu1 = mu[[1]];
		x1 = x[[1]];
		PDF[PoissonDistribution[mu0 + mu1], x1]
	 ] /; TrueQ[Length[x] == Length[mu] == 1]

(* Eq. (37.3), p. 125,  "Discrete Multivariate Distributions" *)
(* bivariate Poisson *)
MultiPoissonDistribution/: PDF[MultiPoissonDistribution[mu0_, mu_?VectorQ,
	aux___], x_?VectorQ] :=
	 Module[{mu1, mu2, x1, x2, i, sum, result},  
		{mu1, mu2} = mu;
		{x1, x2} = x;
		result = Exp[-mu0-mu1-mu2] * sum[(mu1^(x1-i) mu2^(x2-i) mu0^i)/
			( (x1-i)! (x2-i)! i! ), {i, 0, Min[x1, x2]}];
	        result = result /. {sum :> Sum};
		If[!FreeQ[result, Sum],
		    result = result /. {Sum :> sum};
                    result = result /. {i -> Unique[]};
		    result /. {sum :> Sum},
                    result]
	 ] /; TrueQ[Length[x] == Length[mu] == 2]


(* for dimension > 2, use recurrence relations which are useful
	only for numerical x *)
MultiPoissonDistribution/: PDF[MultiPoissonDistribution[mu0_, mu_?VectorQ,
	aux___], x_?VectorQ] :=
	  (
		pP[mu0, mu][x]
	  ) /; TrueQ[Length[x] == Length[mu]]  &&
		If[Apply[And, Map[NumberQ, x]],
		   True,
		   Message[MultiPoissonDistribution::pdf];  False]
		
MultiPoissonDistribution::pdf =
"For dimension greater than two, PDF[MultiPoissonDistribution[mu0, mu], x]
evaluates only for numerical x."

pP[mu0_, mu_][x_ /;
	 !Apply[And, Map[(# == 0)&, x]] && !Apply[Or, Map[(# < 0)&, x]]
	     ] := pP[mu0, mu][x] =
  Module[{i, k = Length[x]},
	 i = First[Position[x, Max[x]]][[1]];
	 mu[[i]]/x[[i]] pP[mu0, mu][x-Table[If[j == i, 1, 0], {j, k}]] +
	 	mu0/x[[i]] pP[mu0, mu][x-Table[1, {k}]]
  ] (* Apply[And, Map[(# >= 1)&, x]] *)

pP[mu0_, mu_][x_] := 0 /; Apply[Or, Map[(# < 0)&, x]]

pP[mu0_, mu_][x_] := (
	Exp[-mu0 -Apply[Plus, mu]]
			) /; Apply[And, Map[(# == 0)&, x]]


(* - - - - - - - - - -  CDF[MultiPoissonDistribution[], ] - - - - - - - - - - *)

(* one or more integers out of range *)
MultiPoissonDistribution/: CDF[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___],
          x_?(VectorQ[#, IntegerQ]&), opts___] := 
	0 /; !Apply[And, Map[(0 <= #)&, x]]


(* general case of CDF *)
MultiPoissonDistribution/:
         CDF[MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___],
          x_?VectorQ, opts___] :=
  Module[{k = Length[mu], y, yvec, prob, arglist, result, flag},
   (
        yvec = Table[Unique[y], {k}];
      (* temporarily turn off the MultiPoissonDistribution::pdf
          message, since we are about to deliberately set up a call
          to the PDF with symbolic args *)
        flag = (Head[MultiPoissonDistribution::pdf] === String);
        Off[MultiPoissonDistribution::pdf];
        prob = PDF[MultiPoissonDistribution[mu0, mu, aux], yvec];
        arglist = Transpose[{yvec, Table[0, {k}], x}];
        PrependTo[arglist, prob];
        result = Apply[Sum, arglist];
        If[flag, On[MultiPoissonDistribution::pdf]];
        If[!FreeQ[result, Sum],
                result /. Thread[Rule[yvec, Table[Unique[], {k}]]],
                result]
   )
  ] /; Apply[And, Map[(!NumberQ[#] || (IntegerQ[#] && 0 <= #))&, x]]


(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* The marginal distributions of the
	 MultiPoissonDistribution[mu0, {mu1, ..., muk}]
	are PoissonDistribution[mu0 + mui]. *)

MultiPoissonDistribution/: Mean[MultiPoissonDistribution[mu0_,
	 mu_?VectorQ, aux___]] :=
	mu0 + mu /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]


(* StandardDeviation yields a vector in multivar case *)
MultiPoissonDistribution/: StandardDeviation[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] :=
   	  Sqrt[mu0 + mu] /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]

(* Variance yields a vector in multivar case *)
MultiPoissonDistribution/: Variance[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] := 
	(mu0 + mu) /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]

MultiPoissonDistribution/: CovarianceMatrix[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] :=
	  With[{k = Length[mu]},
		Table[mu0, {k}, {k}] + DiagonalMatrix[mu]
	  ] /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]

MultiPoissonDistribution/: CorrelationMatrix[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] :=
   Module[{sdvector, sigma},
	  (
		DiagonalMatrix[1/sdvector].sigma.DiagonalMatrix[1/sdvector]
	  ) /; FreeQ[sdvector = StandardDeviation[MultiPoissonDistribution[
		mu0, mu, aux]], StandardDeviation]  &&
	       FreeQ[sigma = CovarianceMatrix[MultiPoissonDistribution[
		mu0, mu, aux]], CovarianceMatrix]	
   ] /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]

MultiPoissonDistribution/: Skewness[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] := 1/Sqrt[mu0 + mu] /;
			ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]

(* A value for the multivariate coefficient of skewness close to 0
	indicates elliptical symmetry. *)
(* From Mardia: A "suitable population counterpart of the measure of
	multivariate skewness for a random vector x distributed
	MultinormalDistribution[mu, sigma] is the expected value of	
	({x - mu}.Inverse[sigma].Transpose[{y-mu}])^3",  where x and
	y are independent and y is also MultinormalDistribution[mu, sigma]).
*)
MultiPoissonDistribution/: MultivariateSkewness[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] :=
   Module[{inverse, k = Length[mu], y, yvec, y0, x, xvec, x0,
	   m = Mean[MultiPoissonDistribution[mu0, mu]],
	   mskewness, table0, rule, table},
    (
	xvec = Table[Unique[x], {k}];
        yvec = Table[Unique[y], {k}];
	(* xvec[[j]] and yvec[[j]] are PoissonDistribution[mu[[j]]],
		j = 1, ..., k,
           and x0 and y0 are each PoissonDistribution[mu0] *)
	mskewness = Expand[
		(({xvec+x0-m}.inverse.Transpose[{yvec+y0-m}])^3)[[1, 1]] ];
	(* Now mskewness is in terms of powers of 2k+2 independent Poisson
		random variables
		x0, y0, xvec[[j]], j = 1, ..., k, and yvec[[j]], j = 1, ..., k.
		Next replace those powers of
		independent random variables with their expectations, i.e.,
		noncentral moments.  The associated Poisson parameter for
		x0 and y0 is mu0.  The associated Poisson parameters for
		xvec and yvec are the elements of mu. *)
	table0 = Table[poissonNoncentralMoment[mu0, j], {j, 3}];
	rule = Join[ Thread[Rule[x0^Range[3], table0]],
		     Thread[Rule[y0^Range[3], table0]] ];
	table = Map[Table[poissonNoncentralMoment[#, j], {j, 3}]&, mu];
	Do[
	   AppendTo[rule, Thread[Rule[ xvec[[i]]^Range[3], table[[i]] ]] ];
	   AppendTo[rule, Thread[Rule[ yvec[[i]]^Range[3], table[[i]] ]] ],
		{i, k}];
	rule = Flatten[rule];
	mskewness /. rule
    ) /; FreeQ[inverse = Inverse[CovarianceMatrix[
		MultiPoissonDistribution[mu0, mu, aux]]], Inverse]
   ] /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]

(* p. 90 of "Discrete Distributions" by Johnson and Kotz *)
poissonNoncentralMoment[theta_, r_] := 
	Sum[StirlingS2[r, j] theta^j, {j, 0, r}]

MultiPoissonDistribution/: Kurtosis[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] :=
		3 + 1/(mu0 + mu) /;
			ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]

(* A value of multivariate kurtosis coefficient close to k*(k+2), where k is
the dimensionality of the data, indicates multinormality. *)
(* From Mardia: A "suitable population counterpart of the measure of
	multivariate kurtosis for a random vector x distributed
	MultinormalDistribution[mu, sigma] is the expected value of	
	({x - mu}.Inverse[sigma].Transpose[{x-mu}])^2
*)
MultiPoissonDistribution/: MultivariateKurtosis[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] := 
   Module[{inverse, k = Length[mu], y, yvec, y0,
	   m = Mean[MultiPoissonDistribution[mu0, mu]],
	   mkurtosis, table0, rule, table },
    (
	yvec = Table[Unique[], {k}];
	(* yvec[[j]] is PoissonDistribution[mu[[j]]], j = 1, ..., k, and
	   y0 is PoissonDistribution[mu0] *)
	mkurtosis = Expand[
	   (({yvec+y0-m}.inverse.Transpose[{yvec+y0-m}])^2)[[1, 1]] ];
	(* Now mkurtosis is in terms of powers of 2k+2 independent Poisson
                random variables
                x0, y0, xvec[[j]], j = 1, ..., k, and yvec[[j]], j = 1, ..., k.
                Next replace those powers of
                independent random variables with their expectations, i.e.,
                noncentral moments.  The associated Poisson parameter for
                x0 and y0 is mu0.  The associated Poisson parameters for
                xvec and yvec are the elements of mu. *)
	table0 = Table[poissonNoncentralMoment[mu0, j], {j, 4}];
        rule = Thread[Rule[y0^Range[4], table0]];
        table = Map[Table[poissonNoncentralMoment[#, j], {j, 4}]&, mu];
        Do[
           AppendTo[rule, Thread[Rule[ yvec[[i]]^Range[4], table[[i]] ]] ],
                {i, k}];
        rule = Flatten[rule];
        mkurtosis /. rule
    ) /; FreeQ[inverse = Inverse[CovarianceMatrix[
		MultiPoissonDistribution[mu0, mu, aux]]], Inverse]
   ] /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]
	
MultiPoissonDistribution/: KurtosisExcess[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] :=
		 1/(mu0 + mu) /;
			ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]

(* A value of multivariate kurtosis excess close to 0 indicates
	multinormality. *) 
MultiPoissonDistribution/: MultivariateKurtosisExcess[MultiPoissonDistribution[
	mu0_, mu_?VectorQ, aux___]] :=
  Module[{k = Length[mu], multikurtosis},
    (
	multikurtosis - k*(k+2)
    ) /; FreeQ[multikurtosis =
		 MultivariateKurtosis[MultiPoissonDistribution[mu0, mu, aux]],
		 MultivariateKurtosis]
  ] /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]
 
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE:
  CharacteristicFunction[dist, t] for a random vector x is defined as the
	expected value of Exp[I t.x], where Length[t] == Length[x].
*)

(* univariate Poisson *)
MultiPoissonDistribution/: CharacteristicFunction[
      MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___], t_?VectorQ] :=
	Module[{mu1, t1},
		mu1 = mu[[1]];
		t1 = t[[1]];	
		CharacteristicFunction[
			PoissonDistribution[mu0 + mu1], t1] 
 	] /; TrueQ[Length[t] == Length[mu] == 1] 
 
(* bivariate Poisson:
	See p. 126 of "Discrete Multivariate Distributions". *)
MultiPoissonDistribution/: CharacteristicFunction[
      MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___], t_?VectorQ] :=
	Module[{mu1, mu2, t1, t2},	
		{mu1, mu2} = mu;
		{t1, t2} = t;
		Exp[-mu1 (1 - Exp[I t1]) - mu2 (1 - Exp[I t2]) -
			mu0 (1 - Exp[I(t1 + t2)]) ]
 	] /; TrueQ[Length[t] == Length[mu] == 2] 
 
(* for dimension > 2, characteristic function is not implemented *)
MultiPoissonDistribution/: CharacteristicFunction[
	MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___], x_?VectorQ] :=
          (
		NotImplemented
          ) /; If[TrueQ[Length[x] == Length[mu] > 2],
		   Message[MultiPoissonDistribution::charfcn];  False,
		   False]

MultiPoissonDistribution::charfcn =
"For dimension greater than two, CharacteristicFunction[
MultiPoissonDistribution[mu0, mu], t] is not implemented."


(*  - - - - - - -  Quantile[MultiPoissonDistribution[], q] - - - - -  - *)
(* NOTE:
Quantile is not implemented for MultiPoissonDistribution[mu0, mu].
Actually since there is no unique vector x such that
q = CDF[MultiPoissonDistribution[mu0, mu], x], it is difficult to even
define Quantile[MultiPoissonDistribution[mu0, mu], q].  One possibility is
the so-called "equicoordinate one-sided percentage points", i.e., the
quantile is {c, c, ..., c} such that CDF[MultiPoissonDistribution[mu0, mu],
{c, c, ..., c}] = q. *)

MultiPoissonDistribution/: Quantile[MultiPoissonDistribution[mu0_, mu_, aux___],
   q_] := NotImplemented /; False

(* - - - -  EllipsoidQuantile[MultiPoissonDistribution[], ellipsoid] - - - - *)
MultiPoissonDistribution/: EllipsoidQuantile[MultiPoissonDistribution[
	mu0_, mu_, aux___], q_] := NotImplemented /; False

(* - - - -  RegionProbability[MultiPoissonDistribution[], ellipsoid] - - - - *)
MultiPoissonDistribution/: RegionProbability[MultiPoissonDistribution[
	mu0_, mu_, aux___], EllipsoidQuantile[mu1_, radii_, dir___?MatrixQ]] :=
			NotImplemented /; False

(* - - - - - - -  ExpectedValue[f, MultiPoissonDistribution[], x] - - - - - - *)
If[FreeQ[Options[ExpectedValue], NSumTerms],
        (* add Options[NSum] *)
   Options[ExpectedValue] = Join[Options[ExpectedValue], Options[NSum]]
]


(* 2 arg syntax *)
MultiPoissonDistribution/: ExpectedValue[f_Function,
  MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___], opts___?OptionQ] :=
  Module[{nn, (* number of arguments of function f *)
	  k = Length[mu], (* number of possible events in a trial; dimension
				of MultiPoissonDistribution[mu0, mu] *)
	  y, yvec, summand, arglist,
          accgoal, compiled, method, extraterms, terms, precgoal,
         verify, wprec, wynn
         },
   (
    yvec = Table[Unique[y], {k}];
    summand = Apply[f, yvec] * PDF[MultiPoissonDistribution[
        mu0, mu, aux], yvec];
    arglist = Map[{#, 0, Infinity}&, yvec];
    PrependTo[arglist, summand];
    {accgoal, compiled, method, extraterms, terms, precgoal,
         verify, wprec, wynn} = {AccuracyGoal, Compiled, Method,
                NSumExtraTerms, NSumTerms, PrecisionGoal,
                VerifyConvergence, WorkingPrecision, WynnDegree} /. {opts} /.
                        Options[ExpectedValue];
    Apply[NSum, Join[arglist,
                {AccuracyGoal->accgoal, Compiled->compiled,
                Method -> method, NSumExtraTerms -> extraterms,
                NSumTerms -> terms, PrecisionGoal -> precgoal,
                VerifyConvergence -> verify, WorkingPrecision -> wprec,
                WynnDegree -> wynn}]
	]
   ) /; (If[Length[f]==1,
	    (* Function with only a body *)
	    nn = Max[Cases[{f}, Slot[z_]->z, Infinity]],
	    (* Function with a list of formal parameters *)
	    nn = Length[f[[1]]]	];
	 nn <= k)
  ] /; ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]


(* 3 arg syntax *)
MultiPoissonDistribution/: ExpectedValue[f_,
	MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___],
		 xvec_?VectorQ, opts___?OptionQ] :=
  Module[{summand, arglist,
          accgoal, compiled, method, extraterms, terms, precgoal,
         verify, wprec, wynn},
    (
	summand = f * PDF[MultiPoissonDistribution[
        	mu0, mu, aux], xvec];
        arglist = Map[{#, 0, Infinity}&, xvec];
        PrependTo[arglist, summand];
        {accgoal, compiled, method, extraterms, terms, precgoal,
         	verify, wprec, wynn} = {AccuracyGoal, Compiled, Method,
                NSumExtraTerms, NSumTerms, PrecisionGoal,
                VerifyConvergence, WorkingPrecision, WynnDegree} /. {opts} /.
                        Options[ExpectedValue];
        Apply[NSum, Join[arglist,
                {AccuracyGoal->accgoal, Compiled->compiled,
                Method -> method, NSumExtraTerms -> extraterms,
                NSumTerms -> terms, PrecisionGoal -> precgoal,
                VerifyConvergence -> verify, WorkingPrecision -> wprec,
                WynnDegree -> wynn}]
	]
    )
  ] /; Length[xvec] <= Length[mu] &&
	ParameterQ[MultiPoissonDistribution[mu0, mu, aux]]
	
	
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE: do not check parameters on Random so that Random will always
        evaluate and be efficient. *)

(* single random vector *)
MultiPoissonDistribution/: Random[MultiPoissonDistribution[mu0_, mu_?VectorQ,
		aux___]] :=
  Module[{ntot, ptot, result},
	Map[Random[PoissonDistribution[#]]&, mu] +
		Random[PoissonDistribution[mu0]]	
  ]	
	     
(* list of random vectors *)
(* NOTE: this RandomArray implementation provides no speed-up tricks such
   as those provided by RandomArray[MultinormalDistribution[]]. *)
MultiPoissonDistribution/: RandomArray[
	MultiPoissonDistribution[mu0_, mu_?VectorQ, aux___], dim_] :=
  Module[{dim1 = dim, arglist, random},
	If[IntegerQ[dim1], dim1 = {dim1}];
	arglist = Prepend[Map[List, dim1],
                        random[MultiPoissonDistribution[mu0, mu, aux]] ];
	Apply[Table, arglist] /. {random :> Random}
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]



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

SetAttributes[ MultinomialDistribution, ReadProtected];
SetAttributes[ NegativeMultinomialDistribution, ReadProtected];
SetAttributes[ MultiPoissonDistribution, ReadProtected];

(* Protect DescriptiveStatistics symbols *)
Protect[Mean, Variance, StandardDeviation, Skewness, Kurtosis,
         KurtosisExcess]

(* Protect MultiDiscreteDistributions symbols. *)
Protect[MultinomialDistribution, NegativeMultinomialDistribution,
	MultiPoissonDistribution]


EndPackage[]


 
