(*:Mathematica Version: 3.0 *)

(*:Package Version: 1.2 *)

(*:Name: Statistics`DiscreteDistributions` *)

(*:Context: Statistics`DiscreteDistributions` *)

(*:Title: Discrete Statistical Distributions *)

(*:Author: David Withoff *)

(*:History:
  David Withoff (Wolfram Research), April 1990.
  Revised March 1991.
  Improved parameter checks for all distributions, added ExpectedValue and
	 RandomArray, speeded up some random number generation,
	 redefined Domain in terms of Range, ECM (Wolfram Research), 1995.
*)

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

(*:Reference: Usage messages only. *)

(*:Summary:
This package provides properties and functionals of discrete probability
distributions used in statistical computations.
*)

(*:Keywords: discrete distribution *)

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

(*:Limitations:
  Currently Sum does not accept the option Assumptions, as does Integrate.
  When Sum is revised to handle assumptions, the ExpectedValue rules 
  should be updated to make the needed assumptions about parameters whereever
  doing so will result in a simplified result.
*)

(*:Warning:
  This package extends the definition of several descriptive
  statistics functions and the definition of Random.  If the original
  usage message for such a function is reloaded, this change will
  not be reflected in the usage message, although the extended
  functionality will remain.
*)

(*:Sources:
  Norman L. Johnson and Samuel Kotz, Discrete Distributions, 1969.
  Stephen Kokoska and Christopher Nevison, Statistical Tables and
    Formulae, 1989.
*)

BeginPackage["Statistics`DiscreteDistributions`",
             "Statistics`DescriptiveStatistics`",
	     "Statistics`Common`DistributionsCommon`"]

(* Extend the usage of Random.  Look for the indicated phrase
        to determine if this has already been done.
*)

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

Random::usage = Random::usage <> " " <>
"Random[distribution] gives a random number with the specified
statistical distribution."

]


 	(* Distributions *)

BernoulliDistribution::usage =
"BernoulliDistribution[p] represents the Bernoulli distribution with
probability p.  A BernoulliDistribution[p] random variable takes on the
values 1 and 0 with probabilities p and 1-p, respectively."

BinomialDistribution::usage =
"BinomialDistribution[n, p] represents the binomial distribution
for n trials with probability p.  A BinomialDistribution[n, p] random
variable describes the number of successes occuring in n trials,
where the probability of success in each trial is p."

DiscreteUniformDistribution::usage =
"DiscreteUniformDistribution[n] represents the discrete uniform
distribution with n possible outcomes. A DiscreteUniformDistribution[n]
random variable takes on the values 1 through n with equal probability." 

GeometricDistribution::usage =
"GeometricDistribution[p] represents the geometric distribution
with probability p.  A GeometricDistribution[p] random variable describes
the number of trials occuring before the first success occurs, where
the probability of success in each trial is p."

HypergeometricDistribution::usage =
"HypergeometricDistribution[n, nsucc, ntot] represents the
hypergeometric distribution for a sample of size n, drawn without
replacement from a population with nsucc successes and total size
ntot.  A HypergeometricDistribution[n, nsucc, ntot] random variable
describes the number of successes occuring in the sample."

LogSeriesDistribution::usage =
"LogSeriesDistribution[theta] represents the logarithmic series
distribution with parameter theta."

NegativeBinomialDistribution::usage =
"NegativeBinomialDistribution[n, p] represents the negative binomial
distribution for success count n and probability p.  A
NegativeBinomialDistribution[n, p] random variable describes the number of
failures that occur before n successes occur, where the probability of
success in each trial is p."

PoissonDistribution::usage =
"PoissonDistribution[mu] represents the Poisson distribution with mean mu."

(*
  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, "specified statistical distribution"] === {},

ExpectedValue::usage = ExpectedValue::usage <> " " <>
"ExpectedValue[f, distribution] gives the expected value of the pure function
f with respect to the specified statistical distribution.
ExpectedValue[f, distribution, x] gives the expected value of the
function f of x with respect to the specified univariate statistical
distribution for x.";
 
Mean::usage = Mean::usage <> " " <>
"Mean[distribution] gives the mean of the specified statistical distribution.";
 
Variance::usage = Variance::usage <> " " <>
"Variance[distribution] gives the variance of the specified statistical
distribution.";
 
StandardDeviation::usage = StandardDeviation::usage <> " " <>
"StandardDeviation[distribution] gives the standard deviation of the
specified statistical distribution.";
 
Skewness::usage = Skewness::usage <> " " <>
"Skewness[distribution] gives the coefficient of skewness of the
specified statistical distribution.";
 
Kurtosis::usage = Kurtosis::usage <> " " <>
"Kurtosis[distribution] gives the coefficient of kurtosis of the
specified statistical distribution.";
 
KurtosisExcess::usage = KurtosisExcess::usage <> " " <>
"KurtosisExcess[distribution] gives the kurtosis excess for the
specified statistical distribution.";
 
Quantile::usage = Quantile::usage <> " " <>
"Quantile[distribution, q] gives the q-th quantile of the specified
statistical distribution."
]
 
Begin["`Private`"]

(* Someday Sum will take the Assumptions option, so allow ExpectedValue to
	take this option as well, even though it is not effective at doing
	anything yet! *)
Options[ExpectedValue] = {Assumptions -> {}}

ExpectedValue::sum =
"In the case of discrete distributions, the Assumptions option of
ExpectedValue is ignored.  Sum does not currently accept Assumptions."


(* ========================================================================= *)
(* Bernoulli Distribution *)

BernoulliDistribution/: ParameterQ[BernoulliDistribution[p_]] := And[
        If[FreeQ[N[p], Complex], True,
           Message[BernoulliDistribution::realparm, p]; False],
        If[0 <= N[p] <= 1, True,
           Message[BernoulliDistribution::probparm, p]; False, True]
]

BernoulliDistribution::probparm =
"Parameter `1` is expected to be a probability between 0 and 1."

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

(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)
BernoulliDistribution/: DomainQ[BernoulliDistribution[p_], list_?VectorQ] :=
	Scan[If[!TrueQ[#==0 || #==1], Return[False]]&, list] =!= False
BernoulliDistribution/: DomainQ[BernoulliDistribution[p_], x_] :=
	TrueQ[x==0 || x==1]

BernoulliDistribution/: Domain[BernoulliDistribution[p_]]:={0, 1} /;
	ParameterQ[BernoulliDistribution[p]]

(* out of range or nonintegral *)
BernoulliDistribution/: PDF[BernoulliDistribution[p_], x_] := 0 /;
        NumberQ[N[x]] && !TrueQ[x==0 || x==1]
(* general case *)
BernoulliDistribution/: PDF[BernoulliDistribution[p_], x_] :=
   If[x == 0, 1-p, p] /;
        ParameterQ[BernoulliDistribution[p]] && Head[x] =!= List

(* zero sum *)
BernoulliDistribution/: CDF[BernoulliDistribution[p_], x_] := 0 /;
        ParameterQ[BernoulliDistribution[p]] && NumberQ[x] && x < 0
(* unit sum *)
BernoulliDistribution/: CDF[BernoulliDistribution[p_], x_] := 1 /;
        ParameterQ[BernoulliDistribution[p]] && NumberQ[x] && x >= 1
(* general case *)
BernoulliDistribution/: CDF[BernoulliDistribution[p_], x_] := 1-p /;
        ParameterQ[BernoulliDistribution[p]] && Head[x] =!= List

BernoulliDistribution/:
    Mean[BernoulliDistribution[p_]] := p /; ParameterQ[
		BernoulliDistribution[p]]

BernoulliDistribution/:
    Variance[BernoulliDistribution[p_]] := p (1-p) /; ParameterQ[
		BernoulliDistribution[p]]

BernoulliDistribution/: StandardDeviation[BernoulliDistribution[p_]] :=
    Sqrt[p (1-p)] /; ParameterQ[BernoulliDistribution[p]]

BernoulliDistribution/: Skewness[BernoulliDistribution[p_]] :=
    (1 - 2 p)/Sqrt[p (1 - p)] /; ParameterQ[BernoulliDistribution[p]]

BernoulliDistribution/: Kurtosis[BernoulliDistribution[p_]] :=
    3 + (1 - 6 p (1-p))/(p (1-p)) /; ParameterQ[BernoulliDistribution[p]]

BernoulliDistribution/: KurtosisExcess[BernoulliDistribution[p_]] :=
    (1 - 6 p (1-p))/(p (1-p)) /; ParameterQ[BernoulliDistribution[p]]

BernoulliDistribution/: CharacteristicFunction[
        BernoulliDistribution[p_], t_] := 1 - p + p Exp[I t] /; ParameterQ[
				BernoulliDistribution[p]]

BernoulliDistribution/: Quantile[BernoulliDistribution[p_], q_] :=
    Which[q <= 1-p, 0, q > 1-p, 1] /; QuantileQ[q] && ParameterQ[
				BernoulliDistribution[p]]

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

bernoulli = Compile[{{p, _Real}}, If[Random[] > p, 0, 1]]

BernoulliDistribution/: Random[BernoulliDistribution[p_]] := bernoulli[p]

BernoulliDistribution/: RandomArray[BernoulliDistribution[p_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[bernoulli[p], {m}];
    If[VectorQ[dim] && Length[dim] > 1,
           Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
           array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

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

(* NOTE: Someday add Assumptions here regarding parameter p *)
BernoulliDistribution/: ExpectedValue[f_Function, BernoulliDistribution[p_],
	opts___?OptionQ] :=
  (If[{opts} =!= {}, Message[ExpectedValue::sum]];
   Sum[f[x] PDF[BernoulliDistribution[p], x], {x, 0, 1}]
  )  /; ParameterQ[BernoulliDistribution[p]]

BernoulliDistribution/: ExpectedValue[f_, BernoulliDistribution[p_], x_Symbol,
	opts___?OptionQ] :=
  (If[{opts} =!= {}, Message[ExpectedValue::sum]];
   Sum[f PDF[BernoulliDistribution[p], x], {x, 0, 1}]
  )  /;		ParameterQ[BernoulliDistribution[p]]


(* ========================================================================= *)
(* Binomial Distribution, K & J, DD, Chap. 3 *)

BinomialDistribution/: ParameterQ[BinomialDistribution[n_, p_]] := And[
    If[FreeQ[N[n], Complex], True,
       Message[BinomialDistribution::realparm, n]; False],
    If[FreeQ[N[p], Complex], True,
       Message[BinomialDistribution::realparm, p]; False],
    If[N[n] > 0, If[IntegerQ[n], True,
        Message[BinomialDistribution::posint, n]; False],
        Message[BinomialDistribution::posint, n]; False, True],
    If[0 <= N[p] <= 1, True,
        Message[BinomialDistribution::probparm, p]; False, True]
]

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

BinomialDistribution::probparm =
"The parameter `1` is expected to be a probability between 0 and 1."

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

(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)
BinomialDistribution/: DomainQ[BinomialDistribution[n_, p_], list_?VectorQ] :=
 Scan[If[!TrueQ[IntegerQ[#] && 0 <= # <= n], Return[False]]&, list] =!= False
BinomialDistribution/: DomainQ[BinomialDistribution[n_, p_], x_] :=
	IntegerQ[x] && TrueQ[0 <= x <= n]

BinomialDistribution/: Domain[BinomialDistribution[n_, p_]]:=
	Module[{domain, range = (Head[Range::range] === $Off)},
	  Off[Range::range];
	  domain = Range[0, n];
	  If[!range, On[Range::range]];
	  domain
	] /; ParameterQ[BinomialDistribution[n, p]]

(* out of range or nonintegral *)
BinomialDistribution/: PDF[BinomialDistribution[n_, p_], x_] := 0 /;
        ParameterQ[BinomialDistribution[n, p]] && NumberQ[N[x]] &&
        (!IntegerQ[x] || (x < 0 || x > n))
(* simplification of special cases *)
BinomialDistribution/: PDF[BinomialDistribution[n_, 0], n_] := 0 /;
        ParameterQ[BinomialDistribution[n, 0]]
BinomialDistribution/: PDF[BinomialDistribution[n_, 1], 0] := 0 /;
        ParameterQ[BinomialDistribution[n, 1]]
(* general case *)
BinomialDistribution/: PDF[BinomialDistribution[n_, p_], x_] :=
   Binomial[n, x] p^x (1-p)^(n-x) /;
        ParameterQ[BinomialDistribution[n, p]] && Head[x] =!= List

(* zero sum *)
BinomialDistribution/: CDF[BinomialDistribution[n_, p_], x_] := 0 /;
        ParameterQ[BinomialDistribution[n, p]] && NumberQ[x] && Negative[x]
(* unit sum *)
BinomialDistribution/: CDF[BinomialDistribution[n_, p_], x_] := 1 /;
        ParameterQ[BinomialDistribution[n, p]] && NumberQ[x] && x >= n
(* simplification of special cases *)
BinomialDistribution/: CDF[BinomialDistribution[n_, p_], n_] := 1 /;
        ParameterQ[BinomialDistribution[n, p]]
BinomialDistribution/: CDF[BinomialDistribution[n_, p_], 0] := (1-p)^n /;
        ParameterQ[BinomialDistribution[n, p]]
BinomialDistribution/: CDF[BinomialDistribution[n_, 1], 0] := 0 /;
        ParameterQ[BinomialDistribution[n, 1]]
(* general case *)
BinomialDistribution/: CDF[BinomialDistribution[n_, p_], x_] :=
   BetaRegularized[1 - p, n-Floor[x], Floor[x]+1] /;
        ParameterQ[BinomialDistribution[n, p]] && Head[x] =!= List

BinomialDistribution/: Mean[BinomialDistribution[n_, p_]] := n p /;
       	ParameterQ[BinomialDistribution[n, p]] 

BinomialDistribution/: Variance[BinomialDistribution[n_, p_]] :=
    n p (1-p) /; ParameterQ[BinomialDistribution[n, p]]

BinomialDistribution/:
        StandardDeviation[BinomialDistribution[n_, p_]] :=
    Sqrt[n p (1-p)] /; ParameterQ[BinomialDistribution[n, p]]

BinomialDistribution/: Skewness[BinomialDistribution[n_, p_]] :=
    (1 - 2 p)/Sqrt[n p (1-p)] /; ParameterQ[BinomialDistribution[n, p]]

BinomialDistribution/: Kurtosis[BinomialDistribution[n_, p_]] :=
    3 + (1 - 6 p (1-p))/(n p (1-p)) /; ParameterQ[BinomialDistribution[n, p]]

BinomialDistribution/: KurtosisExcess[BinomialDistribution[n_, p_]] :=
    (1 - 6 p (1-p))/(n p (1-p)) /; ParameterQ[BinomialDistribution[n, p]]

BinomialDistribution/: CharacteristicFunction[
        BinomialDistribution[n_, p_], t_] :=
    (1 - p + p Exp[I t])^n /; ParameterQ[BinomialDistribution[n, p]]

BinomialDistribution/:
        Quantile[BinomialDistribution[n_, p_], q_] :=
    With[{result = iBinomialQuantile[n, p, q]},
        result /; result =!= Fail] /; 
    NumberQ[n] && NumberQ[N[p]] && ParameterQ[BinomialDistribution[n, p]] &&
								 QuantileQ[q]

iBinomialQuantile[n_, p_, q_] := Module[{low, high, mid},
  If[q==1, Return[n]];
  If[N[(1-p)^n] < q,
    low = 0; high = n;
    While[high - low > 1,
        mid = Floor[(high+low)/2];
        If[N[CDF[BinomialDistribution[n, p], mid]] < q,
            low = mid,
            high = mid,
            high = Fail; Break[]
        ]
    ]; high,
    0,
    Fail
  ]
]

BinomialDistribution/:
        Quantile[BinomialDistribution[n_, p_], q_] :=
		n /; ParameterQ[BinomialDistribution[n, p]] && TrueQ[(q-1)==0]

BinomialDistribution/:
        Quantile[BinomialDistribution[n_, p_], q_] :=
		0 /; ParameterQ[BinomialDistribution[n, p]] && TrueQ[q==0]

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

binomial = Compile[{{n, _Real}, {p, _Real}},
    Apply[Plus, Table[If[Random[] > p, 0, 1], {n}]]	]

BinomialDistribution/: Random[BinomialDistribution[n_, p_]] := binomial[n, p]

BinomialDistribution/: RandomArray[BinomialDistribution[n_, p_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[binomial[n, p], {m}];
    If[VectorQ[dim] && Length[dim] > 1,
           Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
           array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&] 

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

(* NOTE: Someday add Assumptions here regarding parameters n and p *)
BinomialDistribution/: ExpectedValue[f_Function,
		 BinomialDistribution[n_, p_], opts___?OptionQ] :=
   Module[{x, sum},
     (
     If[{opts} =!= {}, Message[ExpectedValue::sum]];
     sum
     ) /; (sum = Sum[ Evaluate[f[x] PDF[BinomialDistribution[n, p], x]],
		         Evaluate[{x, 0, n}] ];
	   FreeQ[sum, Sum])
   ] /; ParameterQ[BinomialDistribution[n, p]]

BinomialDistribution/: ExpectedValue[f_, BinomialDistribution[n_, p_],
		 x_Symbol, opts___?OptionQ] :=
  Module[{sum},
    (
    If[{opts} =!= {}, Message[ExpectedValue::sum]];
    sum
    ) /; (sum = Sum[ Evaluate[f PDF[BinomialDistribution[n, p], x]],
			Evaluate[{x, 0, n}] ];
          FreeQ[sum, Sum])
  ] /; ParameterQ[BinomialDistribution[n, p]]


(* ========================================================================= *)
(* Discrete Uniform Distribution, K & J, DD, Chap. 10, p. 239 *)

ParameterQ[DiscreteUniformDistribution[n_]] := And[
        If[FreeQ[N[n], Complex], True,
           Message[DiscreteUniformDistribution::realparm, n]; False],
        If[N[n] > 0, If[IntegerQ[n], True,
           Message[DiscreteUniformDistribution::posint, n]; False],
           Message[DiscreteUniformDistribution::posint, n]; False, True]
]

DiscreteUniformDistribution::posint =
"The parameter `1` is expected to be a positive integer."

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

(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)
DiscreteUniformDistribution/: DomainQ[DiscreteUniformDistribution[n_],
						list_?VectorQ]:=
 Scan[If[!TrueQ[IntegerQ[#] && TrueQ[1 <= # <= n]], Return[False]]&, list] =!=
						False
DiscreteUniformDistribution/: DomainQ[DiscreteUniformDistribution[n_], x_]:=
	IntegerQ[x] && TrueQ[1 <= x <= n]

DiscreteUniformDistribution/: Domain[DiscreteUniformDistribution[n_]]:=
	Module[{domain, range = (Head[Range::range] === $Off)},
          Off[Range::range];
	  domain = Range[n];
          If[!range, On[Range::range]];
          domain
	] /; ParameterQ[DiscreteUniformDistribution[n]]

(* out of range and nonintegral *)
DiscreteUniformDistribution/: PDF[DiscreteUniformDistribution[n_], x_] := 0 /;
        ParameterQ[DiscreteUniformDistribution[n]] &&
        NumberQ[N[x]] && (!IntegerQ[x] || (x < 1 || x > n))
(* general case *)
DiscreteUniformDistribution/: PDF[DiscreteUniformDistribution[n_], x_] :=
        1/n /;
        ParameterQ[DiscreteUniformDistribution[n]] && Head[x] =!= List

(* zero sum *)
DiscreteUniformDistribution/: CDF[DiscreteUniformDistribution[n_], x_] := 0 /;
        ParameterQ[DiscreteUniformDistribution[n]] &&
                NumberQ[x] && NonPositive[x]
(* unit sum *)
DiscreteUniformDistribution/: CDF[DiscreteUniformDistribution[n_], x_] := 1 /;
        ParameterQ[DiscreteUniformDistribution[n]] &&
        NumberQ[x] && x >= n
(* general case *)
DiscreteUniformDistribution/: CDF[DiscreteUniformDistribution[n_], x_] :=
   Floor[x]/n /;
        ParameterQ[DiscreteUniformDistribution[n]] && Head[x] =!= List

DiscreteUniformDistribution/:
        Mean[DiscreteUniformDistribution[n_]] :=
    (n + 1) / 2 /; ParameterQ[DiscreteUniformDistribution[n]]

DiscreteUniformDistribution/:
        Variance[DiscreteUniformDistribution[n_]] :=
    (n^2 - 1)/12 /; ParameterQ[DiscreteUniformDistribution[n]]

DiscreteUniformDistribution/:
        StandardDeviation[DiscreteUniformDistribution[n_]] :=
    Sqrt[(n^2-1)/12] /; ParameterQ[DiscreteUniformDistribution[n]]

DiscreteUniformDistribution/: Skewness[
        DiscreteUniformDistribution[n_]] := 0 /; ParameterQ[
					DiscreteUniformDistribution[n]]

DiscreteUniformDistribution/:
        Kurtosis[DiscreteUniformDistribution[n_]] :=
    3/5 (3 - 4/(n^2-1)) /; ParameterQ[DiscreteUniformDistribution[n]]

DiscreteUniformDistribution/: KurtosisExcess[
        DiscreteUniformDistribution[n_]] :=
    -(12/(n^2-1) + 6)/5 /; ParameterQ[DiscreteUniformDistribution[n]]

DiscreteUniformDistribution/: CharacteristicFunction[
        DiscreteUniformDistribution[n_], t_] :=
  Module[{k, sum},
    If[FreeQ[sum = Sum[Exp[I t k],{k,1,n}]/n, Sum],
       sum,
       sum /. k -> Unique[]
    ]
  ] /; ParameterQ[DiscreteUniformDistribution[n]]

DiscreteUniformDistribution/: Quantile[
        DiscreteUniformDistribution[n_], q_] :=
    With[{nq = n q},
	     If[NumberQ[nq],
		Max[Ceiling[nq], 1],
		Max[Ceiling[N[nq]], 1]
	     ]
    ] /; ParameterQ[DiscreteUniformDistribution[n]] && QuantileQ[q]

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

DiscreteUniformDistribution/: Random[DiscreteUniformDistribution[n_]] :=
    Random[Integer, {1, n}]

DiscreteUniformDistribution/: RandomArray[
	DiscreteUniformDistribution[n_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[Random[Integer, {1, n}], {m}];
    If[VectorQ[dim] && Length[dim] > 1,
           Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
           array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

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

(* NOTE: someday add Assumptions here regarding parameter n *)
DiscreteUniformDistribution/: ExpectedValue[f_Function,
		DiscreteUniformDistribution[n_], opts___?OptionQ] :=
   Module[{x, sum},
    (
     If[{opts}=!={}, Message[ExpectedValue::sum]];
     sum
    ) /; (sum = Sum[ Evaluate[f[x] PDF[DiscreteUniformDistribution[n], x]],
	 		 Evaluate[{x, 1, n}] ];
	  FreeQ[sum, Sum])
   ] /; ParameterQ[DiscreteUniformDistribution[n]]

DiscreteUniformDistribution/: ExpectedValue[f_,
		DiscreteUniformDistribution[n_], x_Symbol, opts___?OptionQ] :=
  Module[{sum},
  (
	If[{opts}=!={}, Message[ExpectedValue::sum]];
	sum
  ) /; (sum = Sum[ Evaluate[f PDF[DiscreteUniformDistribution[n], x]],
	     Evaluate[{x, 1, n}] ];
        FreeQ[sum, Sum])
  ] /; ParameterQ[DiscreteUniformDistribution[n]]


(* ========================================================================= *)
(* Geometric Distribution, K & J, DD, Chap. 5, p. 122 *)

GeometricDistribution/: ParameterQ[GeometricDistribution[p_]] := And[
        If[FreeQ[N[p], Complex], True,
           Message[GeometricDistribution::realparm, p]; False],
        If[0 <= N[p] <= 1, True,
           Message[GeometricDistribution::probparm, p]; False, True]
]

GeometricDistribution::probparm =
"The parameter `1` is expected to be a probability between 0 and 1."

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

(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)
GeometricDistribution/: DomainQ[GeometricDistribution[p_], list_?VectorQ] :=
	Scan[If[!TrueQ[IntegerQ[#] && # >= 0], Return[False]]&, list] =!= False
GeometricDistribution/: DomainQ[GeometricDistribution[p_], x_] :=
	IntegerQ[x] && x >= 0

GeometricDistribution/: Domain[GeometricDistribution[p_]]:=
        Module[{domain, range = (Head[Range::range] === $Off)},
          Off[Range::range];
          domain = Range[0, Infinity];
          If[!range, On[Range::range]];
          domain
        ] /; ParameterQ[GeometricDistribution[p]]

(* out of range and nonintegral *)
GeometricDistribution/: PDF[GeometricDistribution[p_], x_] := 0 /;
        ParameterQ[GeometricDistribution[p]] &&
        (NumberQ[N[x]] && (!IntegerQ[x] || Negative[x]))
(* general case *)
GeometricDistribution/: PDF[GeometricDistribution[p_], x_] :=
   p (1-p)^x /;
        ParameterQ[GeometricDistribution[p]] && Head[x] =!= List

(* zero sum *)
GeometricDistribution/: CDF[GeometricDistribution[p_], x_] := 0 /;
        ParameterQ[GeometricDistribution[p]] &&
                NumberQ[x] && Negative[x]
(* unit sum *)
GeometricDistribution/: CDF[GeometricDistribution[p_], Infinity] := 1 /;
        ParameterQ[GeometricDistribution[p]]
(* general case *)
GeometricDistribution/: CDF[GeometricDistribution[p_], x_] :=
   1 - (1-p)^Floor[x + 1] /;
        ParameterQ[GeometricDistribution[p]] && Head[x] =!= List

GeometricDistribution/:
    Mean[GeometricDistribution[p_]] := 1/p - 1 /; ParameterQ[
					GeometricDistribution[p]]

GeometricDistribution/: Variance[GeometricDistribution[p_]] :=
    (1 - p) / p^2 /; ParameterQ[GeometricDistribution[p]]

GeometricDistribution/: StandardDeviation[GeometricDistribution[p_]] :=
    Sqrt[1 - p] / p /; ParameterQ[GeometricDistribution[p]]

GeometricDistribution/: Skewness[GeometricDistribution[p_]] :=
    (2 - p)/Sqrt[1 - p] /; ParameterQ[GeometricDistribution[p]]

GeometricDistribution/: Kurtosis[GeometricDistribution[p_]] :=
    (p^2 + 6 - 6 p) / (1 - p) + 3 /; ParameterQ[GeometricDistribution[p]]

GeometricDistribution/: KurtosisExcess[GeometricDistribution[p_]] :=
    (p^2 + 6 - 6 p) / (1 - p) /; ParameterQ[GeometricDistribution[p]]

GeometricDistribution/: CharacteristicFunction[
        GeometricDistribution[p_], t_] :=
    p / (1 - (1-p) Exp[I t]) /; ParameterQ[GeometricDistribution[p]]

GeometricDistribution/: Quantile[GeometricDistribution[p_], q_] :=
    With[{result = iGeometricQuantile[p, q]},
        result /; result =!= Fail
    ] /; ParameterQ[GeometricDistribution[p]] && QuantileQ[q]

iGeometricQuantile[p_, q_] :=
   (If[q == 1, Return[Infinity]];
    With[{l = Log[1-q]/Log[1-p]},
      If[NumberQ[N[l, Precision[l]]],
	If[(NumberQ[p] && p < q) || (!NumberQ[p] && N[p] < q),
		Ceiling[N[l, Precision[l]]] - 1, 0, Fail],
	If[N[p < q], Ceiling[N[l]] - 1, 0, Fail]
      ]
    ])

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

geometric = Compile[{{p, _Real}}, Module[{i=0}, While[Random[] > p, i++];  i]]

GeometricDistribution/: Random[GeometricDistribution[p_]] := geometric[p]

GeometricDistribution/: RandomArray[GeometricDistribution[p_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[geometric[p], {m}];
    If[VectorQ[dim] && Length[dim] > 1,
           Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
           array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]
 
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)

(* NOTE: Someday add Assumptions here regarding parameter p *)
GeometricDistribution/: ExpectedValue[f_Function,
		 GeometricDistribution[p_], opts___?OptionQ] :=
   Module[{x, sum},
    (
     If[{opts}=!={}, Message[ExpectedValue::sum]];
     sum
    ) /; (sum = Sum[ Evaluate[f[x] PDF[GeometricDistribution[p], x]],
	 		 Evaluate[{x, 0, Infinity}] ];
	  FreeQ[sum, Sum])
   ] /; ParameterQ[GeometricDistribution[p]]

GeometricDistribution/: ExpectedValue[f_, GeometricDistribution[p_],
		 x_Symbol, opts___?OptionQ] :=
  Module[{sum},
  (
   If[{opts}=!={}, Message[ExpectedValue::sum]];
   sum
  ) /; (sum = Sum[ Evaluate[f PDF[GeometricDistribution[p], x]],
	     Evaluate[{x, 0, Infinity}] ];
        FreeQ[sum, Sum])
  ] /; ParameterQ[GeometricDistribution[p]]


(* ========================================================================= *)
(* Hypergeometric Distribution, K & J, DD, Chap. 6 *)

HypergeometricDistribution/: 
ParameterQ[HypergeometricDistribution[n_, nsucc_, ntot_]] := And[
If[FreeQ[N[n], Complex], True,
   Message[HypergeometricDistribution::realparm, n]; False],
If[FreeQ[N[nsucc], Complex], True,
   Message[HypergeometricDistribution::realparm, nsucc]; False],
If[FreeQ[N[ntot], Complex], True,
   Message[HypergeometricDistribution::realparm, ntot]; False],
If[N[n] >= 0,
    If[IntegerQ[n], True,
       Message[HypergeometricDistribution::nonnegint, n]; False],
    Message[HypergeometricDistribution::nonnegint, n]; False,
    True],
If[N[nsucc] >= 0,
    If[IntegerQ[nsucc], True,
        Message[HypergeometricDistribution::nonnegint, nsucc]; False],
    Message[HypergeometricDistribution::nonnegint, nsucc]; False,
    True],
If[N[ntot] > 0,
    If[IntegerQ[ntot], True,
        Message[HypergeometricDistribution::posint, ntot]; False],
    Message[HypergeometricDistribution::posint, ntot]; False,
    True],
If[N[n] <= N[ntot], True,
    Message[HypergeometricDistribution::ntoobig, n, ntot]; False, True],
If[N[nsucc] <= N[ntot], True,
    Message[HypergeometricDistribution::nsucctoobig, nsucc, ntot]; False, True]
]

HypergeometricDistribution::nonnegint =
"The parameter `1` is expected to be a nonnegative integer."

HypergeometricDistribution::posint =
"The parameter `1` describing the size of the total population is expected to
be a positive integer."

HypergeometricDistribution::ntoobig =
"The number of trials `1` without replacement must be less than
the total number `2`."

HypergeometricDistribution::nsucctoobig =
"The number of possible successes `1` must be less than the total
number `2`."

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

(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)
HypergeometricDistribution/: DomainQ[
	HypergeometricDistribution[n_, nsucc_, ntot_], list_?VectorQ] :=
  Scan[If[!TrueQ[
	IntegerQ[#] && TrueQ[Max[0,n-ntot+nsucc] <= # <= Min[n,nsucc]]
		 ], Return[False]]&, list] =!= False
HypergeometricDistribution/: DomainQ[
	HypergeometricDistribution[n_, nsucc_, ntot_], x_] :=
	IntegerQ[x] && TrueQ[Max[0,n-ntot+nsucc] <= x <= Min[n,nsucc]]

HypergeometricDistribution/:
   Domain[HypergeometricDistribution[n_, nsucc_, ntot_]]:=
	Module[{domain, range = (Head[Range::range] === $Off)},
          Off[Range::range];
          domain = Range[Max[0,n-ntot+nsucc], Min[n,nsucc]];
          If[!range, On[Range::range]];
          domain
        ] /; ParameterQ[HypergeometricDistribution[n,nsucc,ntot]] &&
	!(MatchQ[{n,nsucc,ntot}, {ntot,ntot,ntot}] ||
	  MatchQ[{n,nsucc,ntot}, {ntot,0,ntot}] ||
	  MatchQ[{n,nsucc,ntot}, {0,ntot,ntot}] ||
	  MatchQ[{n,nsucc,ntot}, {0,0,ntot}])

HypergeometricDistribution/: Domain[HypergeometricDistribution[n_, n_, n_]]:=
	{n} /;	ParameterQ[HypergeometricDistribution[n,n,n]]

HypergeometricDistribution/: Domain[HypergeometricDistribution[l_, m_, n_]]:=
	{0} /;	ParameterQ[HypergeometricDistribution[l,m,n]] && 
	(MatchQ[{l,m}, {n,0}] || MatchQ[{l,m}, {0,n}] || MatchQ[{l,m}, {0,0}])

(* out of range and nonintegral *)
HypergeometricDistribution/: PDF[
   HypergeometricDistribution[n_, nsucc_, ntot_], x_] := 0 /;
        ParameterQ[HypergeometricDistribution[n, nsucc, ntot]] &&
        (NumberQ[N[x]] &&
         (!IntegerQ[x] || (x < Max[0,n-ntot+nsucc] || x > Min[n,nsucc])))
(* special cases *)
HypergeometricDistribution/: PDF[
   HypergeometricDistribution[n_, n_, n_], 0] := 0 /;
        ParameterQ[HypergeometricDistribution[n, n, n]]
HypergeometricDistribution/: PDF[
   HypergeometricDistribution[l_, m_, n_], 0] := 1 /;
        ParameterQ[HypergeometricDistribution[l, m, n]] &&
        (MatchQ[{l,m}, {n,0}] || MatchQ[{l,m}, {0,n}] || MatchQ[{l,m}, {0,0}])
HypergeometricDistribution/: PDF[
   HypergeometricDistribution[l_, m_, n_], n_] := 0 /;
        ParameterQ[HypergeometricDistribution[l, m, n]] &&
        (MatchQ[{l,m}, {n,0}] || MatchQ[{l,m}, {0,n}] || MatchQ[{l,m}, {0,0}])
(* general case *)
HypergeometricDistribution/: PDF[
   HypergeometricDistribution[n_, nsucc_, ntot_], x_] :=
           Binomial[nsucc, x] Binomial[ntot-nsucc, n-x]/Binomial[ntot, n] /;
                ParameterQ[HypergeometricDistribution[n, nsucc, ntot]] &&
                If[NumberQ[x],
                   IntegerQ[x] && Max[0,n-ntot+nsucc] <= x <= Min[n,nsucc],
                   Head[x] =!= List]

(* zero sum *)
HypergeometricDistribution/:
   CDF[HypergeometricDistribution[n_, nsucc_, ntot_], x_] := 0 /;
        ParameterQ[HypergeometricDistribution[n, nsucc, ntot]] &&
        NumberQ[x] && x < Max[0,n-ntot+nsucc]
(* unit sum *)
HypergeometricDistribution/:
   CDF[HypergeometricDistribution[n_, nsucc_, ntot_], x_] := 1 /;
        ParameterQ[HypergeometricDistribution[n, nsucc, ntot]] &&
        NumberQ[x] && x >= Min[n,nsucc]
(* special cases *)
HypergeometricDistribution/:
   CDF[HypergeometricDistribution[n_, n_, n_], 0] := 0  /;
                ParameterQ[HypergeometricDistribution[n, n, n]]
HypergeometricDistribution/:
   CDF[HypergeometricDistribution[l_, m_, n_], 0] := 1  /;
        ParameterQ[HypergeometricDistribution[l, m, n]] &&
        (MatchQ[{l,m}, {n,0}] || MatchQ[{l,m}, {0,n}] || MatchQ[{l,m}, {0,0}])
HypergeometricDistribution/:
   CDF[HypergeometricDistribution[l_, m_, n_], n_] := 1  /;
        (l===n || l==0) && (m===n || m==0)
(* general case *)
HypergeometricDistribution/:
   CDF[HypergeometricDistribution[n_, nsucc_, ntot_], x_] :=
     Module[{k, sum},
       If[FreeQ[sum = Sum[ Evaluate[Binomial[nsucc, k] *
			        Binomial[ntot-nsucc, n-k] / Binomial[ntot, n]],
          		   Evaluate[{k,Max[0,n-ntot+nsucc],Min[x,n,nsucc]}] ], 
		Sum],
	  sum,
	  sum /. k->Unique[]
       ]
     ] /; ParameterQ[HypergeometricDistribution[n, nsucc, ntot]] &&
          Head[x] =!= List

HypergeometricDistribution/:
        Mean[HypergeometricDistribution[n_, nsucc_, ntot_]] :=
    n nsucc / ntot /; ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]

HypergeometricDistribution/:
        Variance[HypergeometricDistribution[n_, nsucc_, ntot_]] :=
    (ntot-n)/(ntot-1) n nsucc/ntot (1-nsucc/ntot) /;
        ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]

HypergeometricDistribution/: StandardDeviation[
        HypergeometricDistribution[n_, nsucc_, ntot_]] :=
    Sqrt[(ntot-n)/(ntot-1) n nsucc (ntot-nsucc)]/ntot /;
        ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]

HypergeometricDistribution/: Skewness[
        HypergeometricDistribution[n_, nsucc_, ntot_]] :=
    (ntot-2 nsucc) (ntot-2 n) Sqrt[ntot-1]/
        ((ntot-2) Sqrt[n nsucc (ntot-nsucc) (ntot-n)]) /;
            ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]

HypergeometricDistribution/: Kurtosis[
        HypergeometricDistribution[n_, nsucc_, ntot_]] :=
    ntot^2 (ntot-1)/((ntot-2)(ntot-3) n nsucc (ntot-nsucc)(ntot-n))*
      (ntot (ntot+1) - 6 n (ntot-n) + 3 nsucc (ntot-nsucc)/ntot^2 *
        (ntot^2 (n-2) - ntot n^2 + 6 n (ntot-n))) /;
            ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]

HypergeometricDistribution/: KurtosisExcess[
        HypergeometricDistribution[n_, nsucc_, ntot_]] :=
    Kurtosis[HypergeometricDistribution[n, nsucc, ntot]] - 3 /;
        ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]

HypergeometricDistribution/: CharacteristicFunction[
        HypergeometricDistribution[n_, nsucc_, ntot_], t_] :=
	((ntot-n)! (ntot-nsucc)!/ntot!) *
	Hypergeometric2F1[-n, -nsucc, ntot-nsucc-n+1, E^(I t)] /;
		ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]

HypergeometricDistribution/: Quantile[
        HypergeometricDistribution[n_, nsucc_, ntot_], q_] :=
    With[{result = iHypergeometricQuantile[n, nsucc, ntot, q]},
        result /; result =!= Fail
    ] /; ParameterQ[HypergeometricDistribution[n, nsucc, ntot]] && QuantileQ[q]

HypergeometricDistribution/: Quantile[
        HypergeometricDistribution[n_, nsucc_, ntot_], q_] :=
	 Min[n, nsucc] /; ParameterQ[
		HypergeometricDistribution[n, nsucc, ntot]] && TrueQ[(q-1)==0]

HypergeometricDistribution/: Quantile[
  HypergeometricDistribution[n_, nsucc_, ntot_], q_] :=
   Max[0, n + nsucc - ntot] /; ParameterQ[
		HypergeometricDistribution[n, nsucc, ntot]] && TrueQ[q==0]

iHypergeometricQuantile[n_, nsucc_, ntot_, q_] :=
Module[{kk = 0, sum, count = Max[0,n-ntot+nsucc]},
    If[q == 1, Return[Min[nsucc, n]]];
    sum = PDF[HypergeometricDistribution[n, nsucc, ntot], count];
    If[NumberQ[sum],
        While[sum < q,
            count++;
            sum += Binomial[nsucc, count] *
                    Binomial[ntot-nsucc, n-count]/
                        Binomial[ntot, n]
        ],
    (* else not a number *)
        count = Fail
    ];
    count
]

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

HypergeometricDistribution/:
        Random[HypergeometricDistribution[n_, nsucc_, ntot_]] :=
		Quantile[HypergeometricDistribution[n,nsucc,ntot], Random[]] 

HypergeometricDistribution/:
	RandomArray[HypergeometricDistribution[n_, nsucc_, ntot_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[Random[HypergeometricDistribution[n, nsucc, ntot]], {m}];
    If[VectorQ[dim] && Length[dim] > 1,
           Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
           array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

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

(* NOTE: Someday add Assumptions here regarding parameters n, nsucc, and ntot *)
HypergeometricDistribution/: ExpectedValue[f_Function,
	HypergeometricDistribution[n_, nsucc_, ntot_], opts___?OptionQ] :=
   Module[{x, sum},
    (
     If[{opts}=!={}, Message[ExpectedValue::sum]];
     sum
    ) /; (sum = Sum[ Evaluate[f[x] *
		 	  PDF[HypergeometricDistribution[n, nsucc, ntot], x]],
		 	 Evaluate[{x, Max[0,n-ntot+nsucc], Min[n,nsucc]}] ];
	  FreeQ[sum, Sum])
   ] /;	ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]

HypergeometricDistribution/: ExpectedValue[f_,
   HypergeometricDistribution[n_, nsucc_, ntot_], x_Symbol, opts___?OptionQ] :=
  Module[{sum},
   (
    If[{opts}=!={}, Message[ExpectedValue::sum]];
    sum
   ) /; (sum = Sum[
	     Evaluate[f PDF[HypergeometricDistribution[n, nsucc, ntot], x]],
	     Evaluate[{x, Max[0,n-ntot+nsucc], Min[n,nsucc]}] ];
         FreeQ[sum, Sum])
  ] /; ParameterQ[HypergeometricDistribution[n, nsucc, ntot]]


(* ========================================================================= *)
(* Logarithmic Series Distribution, K & J, DD, Chap. 7 *)

LogSeriesDistribution/: ParameterQ[LogSeriesDistribution[theta_]] := And[
    If[FreeQ[N[theta], Complex], True,
       Message[LogSeriesDistribution::realparm, theta]; False],
    If[0 <= N[theta] < 1, True,
        Message[LogSeriesDistribution::fracparm, theta]; False,
        True]
]

LogSeriesDistribution::fracparm =
"Parameter `1` is expected to be less than 1 and non-negative."

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

(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)
LogSeriesDistribution/: DomainQ[LogSeriesDistribution[theta_], list_?VectorQ] :=
	Scan[If[!TrueQ[IntegerQ[#] && # >= 1], Return[False]]&, list] =!= False
LogSeriesDistribution/: DomainQ[LogSeriesDistribution[theta_], x_] :=
	IntegerQ[x] && x >= 1

LogSeriesDistribution/: Domain[LogSeriesDistribution[theta_]]:=
        Module[{domain, range = (Head[Range::range] === $Off)},
          Off[Range::range];
          domain = Range[Infinity];
          If[!range, On[Range::range]];
          domain
        ] /; ParameterQ[LogSeriesDistribution[theta]]

(* out of range and nonintegral *)
LogSeriesDistribution/: PDF[LogSeriesDistribution[theta_], x_] := 0 /;
        ParameterQ[LogSeriesDistribution[theta]] &&
        (NumberQ[N[x]] && (!IntegerQ[x] || NonPositive[x]))
(* general case *)
LogSeriesDistribution/: PDF[LogSeriesDistribution[theta_], x_] :=
   theta^x / (-x Log[1-theta]) /;
        ParameterQ[LogSeriesDistribution[theta]] && Head[x] =!= List

(* zero sum *)
LogSeriesDistribution/: CDF[LogSeriesDistribution[theta_], x_] := 0 /;
        ParameterQ[LogSeriesDistribution[theta]] &&
        NumberQ[x] && NonPositive[x]
(* unit sum *)
LogSeriesDistribution/: CDF[LogSeriesDistribution[theta_], Infinity] := 1 /;
        ParameterQ[LogSeriesDistribution[theta]]
(* general case *)	
LogSeriesDistribution/: CDF[LogSeriesDistribution[theta_], x_] :=
 Module[{k, sum},
   If[FreeQ[sum = Sum[ Evaluate[theta^k / k], Evaluate[{k, 1, x}] ], Sum],
      sum,
      sum /. k->Unique[]
   ] / (-Log[1-theta])
 ] /; ParameterQ[LogSeriesDistribution[theta]] && Head[x] =!= List

LogSeriesDistribution/: Mean[LogSeriesDistribution[theta_]] :=
    theta/(1 - theta) / (-Log[1-theta]) /;
	 ParameterQ[LogSeriesDistribution[theta]]

LogSeriesDistribution/: Variance[LogSeriesDistribution[theta_]] :=
  Module[{alpha=-1/Log[1-theta]},
	alpha theta (1-alpha theta)/(1-theta)^2
  ]	/; ParameterQ[LogSeriesDistribution[theta]]

LogSeriesDistribution/: StandardDeviation[LogSeriesDistribution[theta_]] :=
  Module[{alpha=-1/Log[1-theta]},
	Sqrt[alpha theta (1-alpha theta)]/(1-theta)
  ]	/; ParameterQ[LogSeriesDistribution[theta]]

LogSeriesDistribution/: Skewness[LogSeriesDistribution[theta_]] :=
  Module[{alpha=-1/Log[1-theta]},
    (1 + theta - 3 alpha theta + 2 alpha^2 theta^2)/(1-alpha theta)/
	Sqrt[alpha theta(1-alpha theta)]
  ]	/; ParameterQ[LogSeriesDistribution[theta]]

LogSeriesDistribution/: Kurtosis[LogSeriesDistribution[theta_]] :=
  Module[{alpha=-1/Log[1-theta]},
    (1 + 4 theta + theta^2 - 4 alpha theta (1+theta) + 6 alpha^2 theta^2 -
	3 alpha^3 theta^3)/(alpha theta)/(1-alpha theta)^2
  ] 	/; ParameterQ[LogSeriesDistribution[theta]]
 
LogSeriesDistribution/: KurtosisExcess[LogSeriesDistribution[theta_]] :=
	Kurtosis[LogSeriesDistribution[theta]] - 3 /; ParameterQ[
					LogSeriesDistribution[theta]]
 
LogSeriesDistribution/:
        CharacteristicFunction[LogSeriesDistribution[theta_], t_] :=
	Log[1 - theta E^(t I)] / Log[1-theta]	/; ParameterQ[
                                        LogSeriesDistribution[theta]]
 
LogSeriesDistribution/:
        Quantile[LogSeriesDistribution[theta_], q_] :=
    With[{result = iLogSeriesQuantile[theta, q]},
        result /; result =!= Fail] /;
             ParameterQ[LogSeriesDistribution[theta]] && QuantileQ[q]

iLogSeriesQuantile[theta_, q_] := Module[{high, low, mid},
    If[q == 1, Return[Infinity]];
    If[N[theta/(-Log[1-theta])] < q,
        high = 2;
        While[N[CDF[LogSeriesDistribution[theta], high]] < q, high *= 2];
        low = high/2;
        While[high - low > 1,
            mid = (high+low)/2;
            If[N[CDF[LogSeriesDistribution[theta], mid]] < q,
                low = mid,
                high = mid,
                low = high
            ]
        ]; high,
    (* else *)
        1,
    (* indeterminate *)
        Fail
    ]
]

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

LogSeriesDistribution/: Random[LogSeriesDistribution[theta_]] :=
       Quantile[LogSeriesDistribution[theta],Random[]]

LogSeriesDistribution/: RandomArray[LogSeriesDistribution[theta_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[Random[LogSeriesDistribution[theta]], {m}];
    If[VectorQ[dim] && Length[dim] > 1,
           Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
           array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
 
LogSeriesDistribution/: ExpectedValue[f_Function,
		LogSeriesDistribution[theta_], opts___?OptionQ] :=
   Module[{x, sum},
    (
     If[{opts} =!= {}, Message[ExpectedValue::sum]];
     sum
    ) /; (sum = Sum[ Evaluate[f[x] PDF[LogSeriesDistribution[theta], x]],
			 Evaluate[{x, 1, Infinity}] ];
	  FreeQ[sum, Sum])
   ] /; ParameterQ[LogSeriesDistribution[theta]]

LogSeriesDistribution/: ExpectedValue[f_, LogSeriesDistribution[theta_],
	 x_Symbol, opts___?OptionQ] :=
  Module[{sum},
  (
   If[{opts} =!= {}, Message[ExpectedValue::sum]];
   sum
  ) /; (sum = Sum[ Evaluate[ f PDF[LogSeriesDistribution[theta], x]],
       		Evaluate[ {x, 1, Infinity}]];
        FreeQ[sum, Sum])
  ] /; ParameterQ[LogSeriesDistribution[theta]]


(* ========================================================================= *)
(* Negative Binomial Distribution, K & J, DD, Chap. 5 *)

NegativeBinomialDistribution/: ParameterQ[
	NegativeBinomialDistribution[n_, p_]] := And[
    If[FreeQ[N[n], Complex], True,
       Message[NegativeBinomialDistribution::realparm, n]; False],
    If[FreeQ[N[p], Complex], True,
       Message[NegativeBinomialDistribution::realparm, p]; False],
    If[N[n] > 0,
        If[IntegerQ[n], True,
            Message[NegativeBinomialDistribution::posint, n]; False],
        Message[NegativeBinomialDistribution::posint, n]; False, True],
    If[0 < N[p] <= 1, True,
        Message[NegativeBinomialDistribution::probparm, p]; False, True]
]

NegativeBinomialDistribution::posint =
"The parameter `1` describing the failure count is expected to be a positive
integer."

NegativeBinomialDistribution::probparm =
"The parameter `1` is expected to be a non-zero probability
between 0 and 1."

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

(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)
NegativeBinomialDistribution/: DomainQ[
	NegativeBinomialDistribution[n_, p_], list_?VectorQ] :=
	Scan[If[!TrueQ[IntegerQ[#] && # >= 0], Return[False]]&, list] =!= False
NegativeBinomialDistribution/: DomainQ[
	NegativeBinomialDistribution[n_, p_], x_] := IntegerQ[x] && x >= 0

NegativeBinomialDistribution/:
  Domain[NegativeBinomialDistribution[n_, p_]]:=
        Module[{domain, range = (Head[Range::range] === $Off)},
          Off[Range::range];
          domain = Range[0, Infinity];
          If[!range, On[Range::range]];
          domain
        ] /; ParameterQ[NegativeBinomialDistribution[n, p]]

(* out of range and nonintegral *)
NegativeBinomialDistribution/: PDF[NegativeBinomialDistribution[n_, p_], x_] :=
   0 /; ParameterQ[NegativeBinomialDistribution[n, p]] &&
        (NumberQ[N[x]] && (!IntegerQ[x] || Negative[x]))
(* general case *)
NegativeBinomialDistribution/: PDF[NegativeBinomialDistribution[n_, p_], x_] :=
   Binomial[x+n-1, n-1] p^n (1-p)^x /;
        ParameterQ[NegativeBinomialDistribution[n, p]] && Head[x] =!= List

(* zero sum *)
NegativeBinomialDistribution/: CDF[NegativeBinomialDistribution[n_, p_], x_] :=
   0 /; ParameterQ[NegativeBinomialDistribution[n, p]] &&
                NumberQ[x] && Negative[x]
(* unit sum *)
NegativeBinomialDistribution/: CDF[NegativeBinomialDistribution[n_, p_],
        Infinity] := 1 /; ParameterQ[NegativeBinomialDistribution[n, p]]
(* general case *)
NegativeBinomialDistribution/: CDF[NegativeBinomialDistribution[n_, p_], x_] :=
   BetaRegularized[p, n, Floor[x] + 1] /;
        ParameterQ[NegativeBinomialDistribution[n, p]] && Head[x] =!= List

NegativeBinomialDistribution/:
        Mean[NegativeBinomialDistribution[n_, p_]] :=
    n (1-p)/p /; ParameterQ[NegativeBinomialDistribution[n, p]]

NegativeBinomialDistribution/: Variance[
                NegativeBinomialDistribution[n_, p_]] :=
        n (1-p) / p^2 /; ParameterQ[NegativeBinomialDistribution[n, p]]

NegativeBinomialDistribution/: StandardDeviation[
                NegativeBinomialDistribution[n_, p_]] :=
        Sqrt[n (1-p)] / p /; ParameterQ[NegativeBinomialDistribution[n, p]]

NegativeBinomialDistribution/: Skewness[
                NegativeBinomialDistribution[n_, p_]] :=
        (2 - p) / Sqrt[n (1-p)] /; ParameterQ[
					NegativeBinomialDistribution[n, p]]

NegativeBinomialDistribution/:
        Kurtosis[NegativeBinomialDistribution[n_, p_]] :=
    3 + (p^2 + 6 - 6 p) / (n (1-p)) /; ParameterQ[
					NegativeBinomialDistribution[n, p]]

NegativeBinomialDistribution/:
        KurtosisExcess[NegativeBinomialDistribution[n_, p_]] :=
    (p^2 + 6 - 6 p) / (n (1-p)) /; ParameterQ[
					NegativeBinomialDistribution[n, p]]

NegativeBinomialDistribution/: CharacteristicFunction[
        NegativeBinomialDistribution[n_, p_], t_] :=
    (p/(1 - (1-p) Exp[I t]))^n /; ParameterQ[NegativeBinomialDistribution[n, p]]

NegativeBinomialDistribution/:
        Quantile[NegativeBinomialDistribution[n_, p_], q_] :=
    With[{result = iNegativeBinomialQuantile[n, p, q]},
        result /; result =!= Fail] /;
            ParameterQ[NegativeBinomialDistribution[n, p]] && QuantileQ[q]

iNegativeBinomialQuantile[n_, p_, q_] := Module[{high, low, mid},
    If[q == 1, Return[Infinity]];
    If[N[p^n] < q,
        high = 1;
        While[N[CDF[NegativeBinomialDistribution[n, p], high]] < q, high *= 2];
        low = high/2;
        While[high - low > 1,
            mid = (high+low)/2;
            If[N[CDF[NegativeBinomialDistribution[n, p], mid]] < q,
                low = mid,
                high = mid,
                low = high
            ]
        ]; high,
    (* else *)
        0,
    (* indeterminate *)
        Fail
    ]
]

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

negativebinomial = Compile[{{r, _Real}, {p, _Real}}, 
	Module[{nsucc = 0, nfail = 0},
	  If[Random[] > p, nfail++, nsucc++];
	  While[nsucc < r, If[Random[] > p, nfail++, nsucc++]];
	  nfail
	] ]

NegativeBinomialDistribution/: Random[
                NegativeBinomialDistribution[r_, p_]] := negativebinomial[r, p]

NegativeBinomialDistribution/: RandomArray[
		NegativeBinomialDistribution[r_, p_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[negativebinomial[r, p], {m}];
    If[VectorQ[dim] && Length[dim] > 1,
           Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
           array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

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

NegativeBinomialDistribution/: ExpectedValue[f_Function,
	NegativeBinomialDistribution[r_, p_], opts___?OptionQ] :=
  Module[{x, sum},
   (
    If[{opts} =!= {}, Message[ExpectedValue::sum]];
    sum
   ) /; (sum =
	   Sum[ Evaluate[f[x] PDF[NegativeBinomialDistribution[r, p], x]],
	      Evaluate[{x, 0, Infinity}]];
	 FreeQ[sum, Sum])
  ] /; ParameterQ[NegativeBinomialDistribution[n, p]]

NegativeBinomialDistribution/: ExpectedValue[f_,
	 NegativeBinomialDistribution[r_, p_], x_Symbol, opts___?OptionQ] :=
  Module[{sum},
   (
    If[{opts} =!= {}, Message[ExpectedValue::sum]];
    sum
   ) /; (sum = Sum[ Evaluate[f PDF[NegativeBinomialDistribution[r, p], x]],
       		Evaluate[{x, 0, Infinity}] ];
         FreeQ[sum, Sum])
  ] /; ParameterQ[NegativeBinomialDistribution[n, p]]


(* ========================================================================= *)
(* Poisson Distribution, K & J, DD, Chap. 4 *)

PoissonDistribution/: ParameterQ[PoissonDistribution[mu_]] := And[
         If[FreeQ[N[mu], Complex], True,
            Message[PoissonDistribution::realparm, mu]; False],
         If[N[mu] > 0, True,
            Message[PoissonDistribution::posparm, mu]; False, True]
]

PoissonDistribution::posparm =
"Parameter `1` is expected to be greater than zero."

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

(* NOTE: do not condition DomainQ on ParameterQ so that DomainQ will always
        evaluate and be efficient. *)
PoissonDistribution/: DomainQ[PoissonDistribution[mu_], list_?VectorQ] := 
	Scan[If[!TrueQ[IntegerQ[#] && # >= 0], Return[False]]&, list] =!= False
PoissonDistribution/: DomainQ[PoissonDistribution[mu_], x_] := 
	IntegerQ[x] && x >= 0

PoissonDistribution/: Domain[PoissonDistribution[mu_]]:=
        Module[{domain, range = (Head[Range::range] === $Off)},
          Off[Range::range];
          domain = Range[0, Infinity];
          If[!range, On[Range::range]];
          domain
        ] /; ParameterQ[PoissonDistribution[mu]]

(* out of range and nonintegral *)
PoissonDistribution/: PDF[PoissonDistribution[mu_], x_] := 0 /;
        ParameterQ[PoissonDistribution[mu]] &&
        (NumberQ[N[x]] && (!IntegerQ[x] || Negative[x]))
(* general case *)
PoissonDistribution/: PDF[PoissonDistribution[mu_], x_] :=
   Exp[-mu] mu^x / x! /;
        ParameterQ[PoissonDistribution[mu]] && Head[x] =!= List

(* zero sum *)
PoissonDistribution/: CDF[PoissonDistribution[mu_], x_] := 0 /;
        ParameterQ[PoissonDistribution[mu]] && NumberQ[x] && Negative[x]
(* unit sum *)
PoissonDistribution/: CDF[PoissonDistribution[mu_], Infinity] := 1 /;
        ParameterQ[PoissonDistribution[mu]]
(* general case *)
PoissonDistribution/: CDF[PoissonDistribution[mu_], x_] :=
   GammaRegularized[Floor[x] + 1, mu] /;
        ParameterQ[PoissonDistribution[mu]] && Head[x] =!= List

PoissonDistribution/: Mean[PoissonDistribution[mu_]] :=
    mu /; ParameterQ[PoissonDistribution[mu]]

PoissonDistribution/: Variance[PoissonDistribution[mu_]] :=
    mu /; ParameterQ[PoissonDistribution[mu]]

PoissonDistribution/: StandardDeviation[PoissonDistribution[mu_]] :=
    Sqrt[mu] /; ParameterQ[PoissonDistribution[mu]]

PoissonDistribution/: Skewness[PoissonDistribution[mu_]] :=
    1/Sqrt[mu] /; ParameterQ[PoissonDistribution[mu]]

PoissonDistribution/: Kurtosis[PoissonDistribution[mu_]] :=
    3 + 1/mu /; ParameterQ[PoissonDistribution[mu]]

PoissonDistribution/: KurtosisExcess[PoissonDistribution[mu_]] :=
    1/mu /; ParameterQ[PoissonDistribution[mu]]

PoissonDistribution/:
        CharacteristicFunction[PoissonDistribution[mu_], t_] :=
    Exp[mu (Exp[I t] - 1)] /; ParameterQ[PoissonDistribution[mu]]

PoissonDistribution/: Quantile[PoissonDistribution[mu_], q_] :=
    With[{result = iPoissonQuantile[mu, q]},
        result /; result =!= Fail
    ] /; ParameterQ[PoissonDistribution[mu]] && QuantileQ[q]
 
iPoissonQuantile[mu_, q_] := Module[{high, low, mid},
    If[q == 1, Return[Infinity]];
    If[N[Exp[-mu]] < q,
        high = 1;
        While[N[CDF[PoissonDistribution[mu], high]] < q, high *= 2];
        low = high/2;
        While[high - low > 1,
            mid = (high+low)/2;
            If[N[CDF[PoissonDistribution[mu], mid]] < q,
                low = mid,
                high = mid,
                low = high
            ]
        ]; high,
    (* else *)
        0,
    (* indeterminate *)
        Fail
    ]
]

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

PoissonDistribution/: Random[PoissonDistribution[mu_]] :=
       Quantile[PoissonDistribution[mu],Random[]]

PoissonDistribution/: RandomArray[PoissonDistribution[mu_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[Random[PoissonDistribution[mu]], {m}];
    If[VectorQ[dim] && Length[dim] > 1,
           Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
           array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

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

PoissonDistribution/: ExpectedValue[f_Function,
		 PoissonDistribution[mu_], opts___?OptionQ] :=
  Module[{x, sum},
   (
    If[{opts} =!= {}, Message[ExpectedValue::sum]];
    sum
   ) /; (sum = Sum[ Evaluate[f[x] PDF[PoissonDistribution[mu], x]],
		        Evaluate[{x, 0, Infinity}]];
	 FreeQ[sum, Sum])
  ] /; ParameterQ[PoissonDistribution[mu]]

PoissonDistribution/: ExpectedValue[f_, PoissonDistribution[mu_],
	 x_Symbol, opts___?OptionQ] :=
 Module[{sum},
  (
   If[{opts} =!= {}, Message[ExpectedValue::sum]];
   sum
  ) /; (sum = Sum[ Evaluate[f PDF[PoissonDistribution[mu], x]],
       		Evaluate[{x, 0, Infinity}] ];
        FreeQ[sum, Sum])
 ] /; ParameterQ[PoissonDistribution[mu]]


End[]

SetAttributes[ BernoulliDistribution , ReadProtected];
SetAttributes[ BinomialDistribution , ReadProtected];
SetAttributes[ DiscreteUniformDistribution , ReadProtected];
SetAttributes[  GeometricDistribution, ReadProtected];
SetAttributes[ HypergeometricDistribution , ReadProtected];
SetAttributes[  LogSeriesDistribution, ReadProtected];
SetAttributes[ NegativeBinomialDistribution , ReadProtected];
SetAttributes[  PoissonDistribution, ReadProtected];

Protect[BernoulliDistribution, BinomialDistribution,
	DiscreteUniformDistribution, GeometricDistribution,
	HypergeometricDistribution, LogSeriesDistribution,
	NegativeBinomialDistribution, PoissonDistribution];

EndPackage[]

