
(*:Version: Mathematica 2.0 *)

(*:Name: Statistics`DiscreteDistributions *)

(*:Title: Discrete Statistical Distributions *)

(*:Author:
  David Withoff (Wolfram Research), April, 1990.
  Revised March, 1991.
*)

(*:Legal:
  Copyright (c) 1990, 1991, Wolfram Research, Inc.
*)

(*:Reference: Usage messages only. *)

(*:Summary:
  Properties and functionals of discrete probability
  distributions used in statistical computations.
*)

(*:Keywords: discrete 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.
*)

(*: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`"]

 	(* Distributions *)

BernoulliDistribution::usage =
"BernoulliDistribution[p] represents the Bernoulli distribution with mean p."

BinomialDistribution::usage =
"BinomialDistribution[n, p] represents the Binomial distribution
for n trials with probability p."

DiscreteUniformDistribution::usage =
"DiscreteUniformDistribution[n] represents the discrete Uniform
distribution with n possible outcomes."

GeometricDistribution::usage =
"GeometricDistribution[p] represents the Geometric distribution
for probability p."

HypergeometricDistribution::usage =
"HypergeometricDistribution[n, nsucc, ntot] represents the
Hypergeometric distribution for samples of size n drawn without
replacement from a population with nsucc successes and total size
ntot."

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

NegativeBinomialDistribution::usage =
"NegativeBinomialDistribution[n, p] represents the Negative Binomial
distribution for failure count n and probability p."

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

	(*   Functions  *)

(* Introduce usage messages for new functions. *)
 
Domain::usage =
"Domain[distribution] gives the domain of the specified distribution."
 
PDF::usage =
"PDF[distribution, x] gives the probability density function
of the specified statistical distribution evaluated at x."
 
CDF::usage =
"CDF[distribution, x] gives cumulative distribution function of distribution
at the point x, defined as the integral of the probability density of the
the distribution from the lowest value in the domain to x."
 
CharacteristicFunction::usage =
"CharacteristicFunction[distribution, t] gives the characteristic
function of the specified statistical distribution as a function of
the variable t."
 
(*
  Extend usage messages for existing functions.  Look for the
  indicated phrase to determine if this has already been done.
*)

If[StringPosition[Mean::usage, "specified statistical distribution"] === {},
 
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[distribution, q] gives the q-th quantile of the specified
statistical distribution."
]
 
(*
  Extend the usage message 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."
]

Begin["`Private`"]

(* Bernoulli Distribution *)

BernoulliQ[p_] := If[0 <= N[p] <= 1, True,
    Message[BernoulliDistribution::bern, p]; False, True]

BernoulliDistribution::bern =
"Bernoulli distribution parameter `1` is expected to be a probability
between 0 and 1."

BernoulliDistribution/: Domain[BernoulliDistribution[___]]:={0, 1}

BernoulliDistribution/: PDF[BernoulliDistribution[p_], x_] :=
	Which[x == 0, 1-p, x == 1, p, NumberQ[x], 0] /; BernoulliQ[p]

BernoulliDistribution/:
        CDF[BernoulliDistribution[p_], x_] :=
    Which[Negative[x], 0, x >= 1, 1, True,  1-p]  /; BernoulliQ[p]

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

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

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

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

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

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

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

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

BernoulliDistribution/: Random[BernoulliDistribution[p_]] :=
        If[Random[] > p, 0, 1] /; BernoulliQ[p]

(* Binomial Distribution *)

BinomialPQ[n_, p_] := And[
    If[N[n] > 0, If[IntegerQ[n], True,
        Message[BinomialDistribution::trials, n]; False],
        Message[BinomialDistribution::trials, n]; False, True],
    If[0 <= N[p] <= 1, True,
        Message[BinomialDistribution::probability, p]; False, True]
]

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

BinomialDistribution/: Domain[BinomialDistribution[n_, p_]]:= Range[n+1]-1 /;
				IntegerQ[n] && 1<=n<=2 && BinomialPQ[n,p]
					
BinomialDistribution/: Domain[BinomialDistribution[n_, p_]]:=
	{0,1," ..." ,n}		/; BinomialPQ[n,p]

BinomialDistribution/: PDF[BinomialDistribution[n_, p_], x_] :=
    If[0 <= x <= n,
        If[IntegerQ[x], Binomial[n, x] p^x (1-p)^(n-x), 0], 0] /;
        BinomialPQ[n, p]

BinomialDistribution/: CDF[BinomialDistribution[n_, p_], x_] :=
    Which[
        Negative[x], 0,
        x >= n, 1,
        True, BetaRegularized[1 - p, n-Floor[x], Floor[x]+1]
    ] /; BinomialPQ[n, p]
         
BinomialDistribution/: Mean[BinomialDistribution[n_, p_]] := n p /;
        BinomialPQ[n, p]

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

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

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

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

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

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

BinomialDistribution/:
        Quantile[BinomialDistribution[n_, p_], q_] :=
    With[{result = iBinomialQuantile[n, p, q]},
        result /; result =!= Fail] /; 
    NumberQ[n] && NumberQ[N[p]] && BinomialPQ[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/: Random[BinomialDistribution[n_, p_]] :=
    Quantile[BinomialDistribution[n, p], Random[]]

(* Discrete Uniform Distribution *)

DiscreteUniformQ[n_] := If[N[n] > 0, If[IntegerQ[n], True,
        Message[DiscreteUniformDistribution::count, n]; False],
        Message[DiscreteUniformDistribution::count, n]; False, True]
 
DiscreteUniformDistribution::count =
"The parameter `1` is expected to be a positive integer."
 
DiscreteUniformDistribution/: Domain[DiscreteUniformDistribution[n_]]:=
	Range[n] /; IntegerQ[n] && 1<=n<=3 && DiscreteUniformQ[n]

DiscreteUniformDistribution/:
	Domain[DiscreteUniformDistribution[n_]]:={1,2, "..." ,n} /;
						DiscreteUniformQ[n]

DiscreteUniformDistribution/:
        PDF[DiscreteUniformDistribution[n_], x_] :=
    If[1 <= x <= n, If[IntegerQ[x], 1/n, 0], 0] /; DiscreteUniformQ[n]

DiscreteUniformDistribution/: CDF[DiscreteUniformDistribution[n_], x_] :=
    Which[
        Negative[x], 0,
        x >= n, 1,
        True, Max[0, Min[1, Floor[x]/n]]
    ] /;	DiscreteUniformQ[n]
 
DiscreteUniformDistribution/:
        Mean[DiscreteUniformDistribution[n_]] :=
    (n + 1) / 2 /; DiscreteUniformQ[n]

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

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

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

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

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

DiscreteUniformDistribution/: CharacteristicFunction[
        DiscreteUniformDistribution[n_], t_] :=
    Sum[Exp[I t k],{k,1,n}]/n	/; DiscreteUniformQ[n]

DiscreteUniformDistribution/: Quantile[
        DiscreteUniformDistribution[n_], q_] :=
    Max[Ceiling[N[n q]], 1] /; DiscreteUniformQ[n] && QuantileQ[q]

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

(* Geometric Distribution *)

GeometricQ[p_] := If[0 <= N[p] <= 1, True,
        Message[GeometricDistribution::probability, p]; False, True]

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

GeometricDistribution/:
	Domain[GeometricDistribution[___]]:= {0,1,2, "..."}

GeometricDistribution/:
        PDF[GeometricDistribution[p_], x_] :=
    If[!Negative[x], If[IntegerQ[x], p (1-p)^x, 0], 0] /; GeometricQ[p]

GeometricDistribution/: 
        CDF[GeometricDistribution[p_], x_] :=
    Which[
        Negative[x], 0,
        x == Infinity, 1,
        True, 1 - (1-p)^Floor[x + 1]
    ] /;	GeometricQ[p]

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

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

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

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

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

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

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

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

iGeometricQuantile[p_, q_] :=
   (If[q == 1, Return[Infinity]];
    If[N[p] < q,
	Ceiling[N[Log[1-q]/Log[1-p]]] - 1,
	0,
	Fail
    ])

GeometricDistribution/: Random[GeometricDistribution[p_]] :=
    Quantile[GeometricDistribution[p], Random[]]

(* Hypergeometric Distribution *)

HypergeometricPQ[n_, nsucc_, ntot_] := And[
If[N[n] >= 0,
    If[IntegerQ[n], True,
        Message[HypergeometricDistribution::posint, n]; False],
    Message[HypergeometricDistribution::posint, n]; False,
    True],
If[N[nsucc] >= 0,
    If[IntegerQ[nsucc], True,
        Message[HypergeometricDistribution::posint, nsucc]; False],
    Message[HypergeometricDistribution::posint, 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::posint =
"The parameter `1` 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/:
	Domain[HypergeometricDistribution[n_, nsucc_, ntot_]]:=
		{Max[0,n-ntot+nsucc], "...", Min[n,nsucc]} /;
			HypergeometricPQ[n,nsucc,ntot]

HypergeometricDistribution/: PDF[
        HypergeometricDistribution[n_, nsucc_, ntot_], x_] :=
    If[Max[0,n-ntot+nsucc] <= x <= Min[n,nsucc],
        If[IntegerQ[x],
            Binomial[nsucc, x] Binomial[ntot-nsucc, n-x]/
                Binomial[ntot, n], 0],
        0] /; HypergeometricPQ[n, nsucc, ntot]

HypergeometricDistribution/:
        CDF[HypergeometricDistribution[n_, nsucc_, ntot_], x_] :=
    Which[
	x < Max[0,n-ntot+nsucc], 0,
	x >= Min[n,nsucc], 1,
	True, Sum[PDF[HypergeometricDistribution[n,nsucc,ntot],k],
			{k,Max[0,n-ntot+nsucc],Min[x,n,nsucc]}]
    ] /; HypergeometricPQ[n, nsucc, ntot]

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

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

HypergeometricDistribution/: StandardDeviation[
        HypergeometricDistribution[n_, nsucc_, ntot_]] :=
    Sqrt[(ntot-n)/(ntot-1) n nsucc (ntot-nsucc)]/ntot /;
        HypergeometricPQ[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)]) /;
            HypergeometricPQ[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))) /;
            HypergeometricPQ[n, nsucc, ntot]

HypergeometricDistribution/: KurtosisExcess[
        HypergeometricDistribution[n_, nsucc_, ntot_]] :=
    Kurtosis[HypergeometricDistribution[n, nsucc, ntot]] - 3 /;
        HypergeometricPQ[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)] /;
		HypergeometricPQ[n, nsucc, ntot]

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

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
]

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

(* Logarithmic Series Distribution *)

LogSeriesQ[theta_] :=
    If[0 <= N[theta] < 1, True,
        Message[LogSeriesDistribution::outofrange, theta]; False,
        True]

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

LogSeriesDistribution/: Domain[LogSeriesDistribution[___]]:= {1,2,3, "..."}

LogSeriesDistribution/: PDF[LogSeriesDistribution[theta_], x_] :=
    If[Positive[x], If[IntegerQ[x], theta^x / (-x Log[1-theta]), 0], 0] /;
							LogSeriesQ[theta]

LogSeriesDistribution/:
        CDF[LogSeriesDistribution[theta_], x_] :=
    Which[
      x < 1, 0,
      x == Infinity, 1,
      True,
   	Sum[theta^k / k, {k, 1, x}] / (-Log[1-theta])
    ]	/; LogSeriesQ[theta]

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

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

LogSeriesDistribution/: StandardDeviation[LogSeriesDistribution[theta_]] :=
  Module[{alpha=-1/Log[1-theta]},
	Sqrt[alpha theta (1-alpha theta)]/(1-theta)
  ]	/; LogSeriesQ[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)]
  ]	/; LogSeriesQ[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
  ] 	/; LogSeriesQ[theta]
 
LogSeriesDistribution/: KurtosisExcess[LogSeriesDistribution[theta_]] :=
	Kurtosis[LogSeriesDistribution[theta]] - 3 /; LogSeriesQ[theta]
 
LogSeriesDistribution/:
        CharacteristicFunction[LogSeriesDistribution[theta_], t_] :=
	Log[1 - theta E^(t I)] / Log[1-theta]	/; LogSeriesQ[theta]
 
LogSeriesDistribution/:
        Quantile[LogSeriesDistribution[theta_], q_] :=
    With[{result = iLogSeriesQuantile[theta, q]},
        result /; result =!= Fail] /;
             LogSeriesQ[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
    ]
]

LogSeriesDistribution/: Random[LogSeriesDistribution[theta_]] :=
       Quantile[LogSeriesDistribution[theta],Random[]] /; LogSeriesQ[theta]
 
(* Negative Binomial Distribution *)

NegativeBinomialPQ[n_, p_] := And[
    If[N[n] > 0,
        If[IntegerQ[n], True,
            Message[NegativeBinomialDistribution::trials, n]; False],
        Message[NegativeBinomialDistribution::trials, n]; False, True],
    If[0 < N[p] <= 1, True,
        Message[NegativeBinomialDistribution::probability, p]; False, True]
]

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

NegativeBinomialDistribution/:
	Domain[NegativeBinomialDistribution[___]]:={0,1,2, "..."}

NegativeBinomialDistribution/:
        PDF[NegativeBinomialDistribution[n_, p_], x_] :=
    If[!Negative[x],
        If[IntegerQ[x], Binomial[x+n-1, n-1] p^n (1-p)^x, 0],
        0] /; NegativeBinomialPQ[n, p]

NegativeBinomialDistribution/:
        CDF[NegativeBinomialDistribution[n_, p_], x_] :=
    Which[
      Negative[x], 0,
      x == Infinity, 1,
      True,
	BetaRegularized[p, n, Floor[x] + 1]
    ]	/; NegativeBinomialPQ[n, p]

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

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

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

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

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

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

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

NegativeBinomialDistribution/:
        Quantile[NegativeBinomialDistribution[n_, p_], q_] :=
    With[{result = iNegativeBinomialQuantile[n, p, q]},
        result /; result =!= Fail] /;
            NegativeBinomialPQ[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
    ]
]

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

(* Poisson Distribution *)

PoissonPQ[mu_] := If[Positive[mu], True,
    Message[PoissonDistribution::parameter, mu]; False, True]

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

PoissonDistribution/: Domain[PoissonDistribution[___]]:={0,1,2, "..."}

PoissonDistribution/: PDF[PoissonDistribution[mu_], x_] :=
    If[!Negative[x], If[IntegerQ[x], Exp[-mu] mu^x / x!, 0], 0] /;
        PoissonPQ[mu]

PoissonDistribution/: CDF[PoissonDistribution[mu_], x_] :=
    Which[
      Negative[x], 0,
      x == Infinity, 1,
      True,
	GammaRegularized[Floor[x] + 1, mu, Infinity]
    ]	/; PoissonPQ[mu]

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

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

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

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

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

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

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

PoissonDistribution/: Quantile[PoissonDistribution[mu_], q_] :=
    With[{result = iPoissonQuantile[mu, q]},
        result /; result =!= Fail
    ] /; PoissonPQ[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
    ]
]

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

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[]

