(*:Mathematica Version: 4.0 *)

(*:Package Version: 1.3 *)

(*:Name: Statistics`ContinuousDistributions` *)

(*:Context: Statistics`ContinuousDistributions` *)

(*:Title: Continuous Statistical Distributions *)

(*:Author: David Withoff *)

(*:History:
  David Withoff (Wolfram Research), February 1990.
  Modified March 1991.
  Improved efficiency of CDF of noncentral distributions,
	Jeff Adams (Wolfram Research), February 1993.
  Removed import of Statistics`InverseStatisticalFunctions` because the inverse
	  functions of Erf, Erfc, GammaRegularized, and BetaRegularized
	  are autoloaded now, ECM (Wolfram Research), October 1994.
  Added ParameterQ, DomainQ, ExpectedValue, RandomArray, and
	ParetoDistribution, speeded up random number generation,
	redefined Domain in terms of Interval,
	ECM (Wolfram Research), February 1995.
*)

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

(*:Reference: Usage messages only. *)

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

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

  Definitions for ChiSquareDistribution, FRatioDistribution,
  NormalDistribution, and StudentTDistribution are contained
  in the Statistics`NormalDistribution` package.
*)

(*:Limitations:
CharacteristicFunction is not implemented for
  NoncentralChiSquareDistribution, NoncentralFRatioDistribution,
  NoncentralStudentTDistribution, LogNormalDistribution, and 
  WeibullDistribution.
*)

(*:Note: Most functions provide numerical results.  The CDF of noncentral
	distributions NoncentralChiSquareDistribution,
	NoncentralFRatioDistribution, and NoncentralStudentTDistribution,
	is calculated via NIntegrate.  To calculate the CDF with a higher
	working precision, use SetOptions[NIntegrate, WorkingPrecision->val]
	where val is greater than the default value given by
	Options[NIntegrate].
*)
 
(*:Sources:
  Norman L. Johnson and Samuel Kotz, Continuous Univariate
    Distributions, 1970, Volumes 1 & 2, John Wiley & Sons.
  M.E.Johnson, Multivariate Statistical Simulation, 1987, John Wiley & Sons.
  R.C.H.Cheng & G.M.Feast, "Some Simple Gamma Variate Generators", 
   Applied Statistics (1979), 28, No. 3, pp. 290-295. 
*)


BeginPackage["Statistics`ContinuousDistributions`",
             "Statistics`NormalDistribution`",
             "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."

]

(* Quantile needs to accept options for
   NoncentralChiSquareDistribution, NoncentralFRatioDistribution and
   NoncentralStudentTDistribution.  Add those options if needed. *)

If[FreeQ[Options[Quantile], WorkingPrecision],

Unprotect[Quantile];
Options[Quantile] = Join[Options[Quantile], 
	{WorkingPrecision -> $MachinePrecision,
	 AccuracyGoal -> Automatic, PrecisionGoal -> Automatic,
	 MaxIterations -> 30}];
Protect[Quantile]

]

	(*  Distributions  *)

BetaDistribution::usage =
"BetaDistribution[p, q] represents the continuous beta distribution
with shape parameters p and q."

CauchyDistribution::usage =
"CauchyDistribution[a, b] represents the Cauchy distribution with
location parameter a and scale parameter b."

ChiDistribution::usage =
"ChiDistribution[n] represents the chi distribution with n degrees
of freedom."
 
NoncentralChiSquareDistribution::usage =
"NoncentralChiSquareDistribution[n, lambda] represents the non-central
chi-square distribution with n degrees of freedom and non-centrality
parameter lambda."
 
ExponentialDistribution::usage =
"ExponentialDistribution[lambda] represents the exponential
distribution with scale inversely proportional to parameter lambda."
 
ExtremeValueDistribution::usage =
"ExtremeValueDistribution[alpha, beta] represents the extreme-value
(Fisher-Tippett) distribution with location parameter alpha and scale
parameter beta."

NoncentralFRatioDistribution::usage =
"NoncentralFRatioDistribution[n1, n2, lambda] represents the non-central
F distribution with n1 numerator degrees of freedom, n2 denominator
degrees of freedom, and numerator non-centrality parameter lambda."

GammaDistribution::usage =
"GammaDistribution[alpha, beta] represents the gamma distribution
with shape parameter alpha and scale parameter beta."

HalfNormalDistribution::usage =
"HalfNormalDistribution[theta] represents the half-normal distribution
with scale inversely proportional to parameter theta."

LaplaceDistribution::usage =
"LaplaceDistribution[mu, beta] represents the Laplace (double
exponential) distribution with mean mu and scale parameter beta."

LogNormalDistribution::usage =
"LogNormalDistribution[mu, sigma] represents the log-normal distribution
based on a normal distribution having mean mu and standard deviation sigma."

LogisticDistribution::usage =
"LogisticDistribution[mu, beta] represents the logistic distribution
with mean mu and scale parameter beta."

RayleighDistribution::usage =
"RayleighDistribution[sigma] represents the Rayleigh distribution
with scale parameter sigma."
 
NoncentralStudentTDistribution::usage =
"NoncentralStudentTDistribution[n, delta] represents the non-central
Student's t distribution with n degrees of freedom and non-centrality
parameter delta."

ParetoDistribution::usage =
"ParetoDistribution[k, alpha] represents the Pareto distribution with 
minimum value parameter k and shape parameter alpha."

UniformDistribution::usage =
"UniformDistribution[min, max] represents the uniform distribution
on the interval {min, max}."

WeibullDistribution::usage =
"WeibullDistribution[alpha, beta] represents the Weibull distribution
with shape parameter alpha and scale parameter beta."

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

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."

]

Unprotect[BetaDistribution,CauchyDistribution,ChiDistribution,
	NoncentralChiSquareDistribution,ExponentialDistribution, 
	ExtremeValueDistribution,NoncentralFRatioDistribution,
	GammaDistribution,HalfNormalDistribution, LaplaceDistribution,
	LogNormalDistribution,LogisticDistribution,
	RayleighDistribution,NoncentralStudentTDistribution, 
	ParetoDistribution, UniformDistribution, WeibullDistribution];

Begin["`Private`"]

quantileOptionCheck[workprec_, accgoal_, precgoal_, maxits_] :=
  Module[{pg = precgoal, ag = accgoal, wp = workprec}, 
      If[!(IntegerQ[maxits] && Positive[maxits]),
           Message[Quantile::badits, maxits];
           maxits = 30];
      If[!(IntegerQ[wp] && wp >= $MachinePrecision),
           Message[Quantile::badwork, wp];
           wp = $MachinePrecision];
      If[pg === Automatic, pg = wp - 10];
      If[ag === Automatic, ag = wp - 10];
      If[!(IntegerQ[pg] && Positive[pg]),
           Message[Quantile::badprec, pg];
           pg = $MachinePrecision - 10];
      If[!(IntegerQ[ag] && Positive[ag]),
           Message[Quantile::badacc, ag];
           ag = $MachinePrecision - 10];
      {wp, ag, pg, maxits}
  ]

(* NOTE: this could be speeded up by using NIntegrate[f, {y, y0, y1}] to
	help calculate NIntegrate[f, {y, y0, y2}], y0 < y1 < y2.  However,
	doing this complicates the tracking of accuracy and precision. *)
findQuantile[pdf_, {x_, x0_, x1_}, {qqmin_, qqmax_},
		 qfraction_, quantile0_, ag_, pg_, wp_, maxits_] :=
  Module[{quantile = quantile0, qmin = qqmin, qmax = qqmax,
	  newquantile, pdfAtQuantile, limitlist, int, delta,
          q = SetPrecision[qfraction, wp],
	  accgoal = N[10^-ag], precgoal = N[10^-pg],
	  accDelta, NAcc, precDelta, NPrec, its = 1,
	  accflag0, accflag1 = False,        		(* flags *)
	  precflag0, precflag1 = False,     		(* flags *)
	  invaliddigits, validrange,                    (* flags *)
	  ovflag = False,				(* flags *) 
	  precw, accbd, accsm, mnprc  (* messages *) },
   pdfAtQuantile = pdf /. {x->quantile};
   limitlist = {x, x0, x1} /. Null -> quantile;
   precw = Head[NIntegrate::precw]===$Off;  Off[NIntegrate::precw];
   int = NIntegrate[Evaluate[pdf], Evaluate[limitlist],
		Evaluate[WorkingPrecision -> wp]];
   If[x0 === Null, If[int < qfraction, qmax = quantile, qmin = quantile]];
   If[x1 === Null, If[int < qfraction, qmin = quantile, qmax = quantile]];
   delta = (int - q)/pdfAtQuantile;
   newquantile = quantile - delta;
   If[newquantile < qmin, newquantile = (quantile + qmin)/2];
   If[newquantile > qmax, newquantile = (quantile + qmax)/2];
   ovflag = Check[
      While[(
	      (* Check that there are enough digits of accuracy in
                 accDelta to estimate accuracy of quantile, and that
                 there are enough digits of precision in precDelta
                 to estimate precision of quantile. *)
              (	
                accDelta = Abs[quantile - newquantile];
	        NAcc = -Floor[Log[10, accDelta]]-1;
                precDelta = accDelta/(Abs[quantile] + Abs[newquantile])/2;
	        NPrec = -Floor[Log[10, precDelta]]-1;
	        invaliddigits = TrueQ[
		 ((NAcc =!= Infinity) && (Accuracy[accDelta] < NAcc)) &&
                 ((NPrec =!= Infinity) && (Precision[precDelta] < NPrec))];
	        !invaliddigits ) &&

	      ( accflag0 = (0 <= accDelta < accgoal);
	        precflag0 = (0 <= precDelta < precgoal); 
	        !((accflag0 && accflag1) || (precflag0 && precflag1)) ) &&

	      its <= maxits &&
	   
	      (validrange = (qqmin <= quantile <= qqmax))
	    ),

        {accflag1, precflag1} = {accflag0, precflag0};

        quantile = newquantile;
        pdfAtQuantile = pdf /. {x->quantile};
        limitlist = {x, x0, x1} /. Null -> quantile;
	int = NIntegrate[Evaluate[pdf], Evaluate[limitlist],
			Evaluate[WorkingPrecision -> wp]];
        If[x0 === Null, If[int < qfraction, qmax = quantile, qmin = quantile]];
        If[x1 === Null, If[int < qfraction, qmin = quantile, qmax = quantile]];
        delta = (int - q)/pdfAtQuantile;
	newquantile = quantile - delta;
        If[newquantile < qmin, newquantile = (quantile + qmin)/2];
        If[newquantile > qmax, newquantile = (quantile + qmax)/2];
        its++

      ], (* end While *)
      True, General::ovfl];   (* end Check *)
      
      If[!precw, On[NIntegrate::precw]];

      quantile = newquantile;
      
      (* check other error possibilities *)
      If[!validrange, Message[Quantile::outofrng];  Return[$Failed]];
      If[its > maxits, Message[Quantile::qcv, maxits]];
      If[invaliddigits,
            If[wp === $MachinePrecision,
               Message[Quantile::qmp],
               Message[Quantile::qdig, wp]]
     ];
     If[TrueQ[ovflag], Message[Quantile::over]];

     (* Adjust the accuracy and precision of quantile. *)
     accbd = Head[SetAccuracy::accbd]===$Off;
     accsm = Head[SetAccuracy::accsm]===$Off;
     mnprc = Head[Precision::mnprc]===$Off;
     Off[SetAccuracy::accbd];  Off[SetAccuracy::accsm];  Off[Precision::mnprc];
     If[Accuracy[quantile] > NAcc,
            quantile = SetAccuracy[quantile, NAcc]	];
     If[Precision[quantile] > NPrec,
            quantile = SetPrecision[quantile, NPrec]	];
     If[!accbd, On[SetAccuracy::accbd]];  If[!accsm, On[SetAccuracy::accsm]];
     If[!mnprc, On[Precision::mnprc]];

     quantile

  ] (* end findQuantile *)

Quantile::over =
"Warning: numeric overflow occurred during this search.  This may mean that
the search is starting from an inappropriate point or that
insufficient precision is being used for these calculations.
The returned parameter quantile may not be correct.";

Quantile::qcv =		(* i.e., similar to FindMinimum::fmcv *)
"Quantile failed to converge to the requested accuracy or precision for
the quantile within `1` iterations.";

Quantile::qdig =	(* i.e., similar to FindMinimum::fmdig *)
"`1` working digits is insufficient to achieve the requested accuracy or
precision for the quantile.";

Quantile::qmp =		(* i.e., similar to FindMinimum::fmmp *)
"Machine precision is insufficient to achieve the requested accuracy or
precision for the sum of squares.";

Quantile::outofrng =
"The search for the quantile exited the allowable range.";

Quantile::badits =
"Warning: MaxIterations -> `1` is not set to a positive integer;
setting to 30.";

Quantile::badwork =
"Warning: The value `1` given for the WorkingPrecision is not an integer or is
an integer less than $MachinePrecision.  Setting to $MachinePrecision.";

Quantile::badprec =
"Warning: The value `1` given for the PrecisionGoal is not a positive integer
or Automatic.  Setting to $MachinePrecision - 10.";

Quantile::badacc =
"Warning: The value `1` given for the AccuracyGoal is not a positive integer
or Automatic.  Setting to $MachinePrecision - 10.";


(* ======================== Beta Distribution ========================= *)
(* Continuous Beta (Binomial) Distribution, J & K, Vol.2, Chap. 24 *)

BetaDistribution/: ParameterQ[BetaDistribution[p_, q_]] := And[
        If[FreeQ[N[p], Complex], True,
           Message[BetaDistribution::realparm, p]; False],
        If[FreeQ[N[q], Complex], True,
           Message[BetaDistribution::realparm, q]; False],
        If[N[p] > 0, True,
           Message[BetaDistribution::posparm, p]; False, True],
        If[N[q] > 0, True,
           Message[BetaDistribution::posparm, q]; False, True]
]

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

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

BetaDistribution/: DomainQ[BetaDistribution[p_, q_], list_?VectorQ] := 
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[0 <= # <= 1], Return[False]]&, list] =!= False
BetaDistribution/: DomainQ[BetaDistribution[p_, q_], x_] := 
	FreeQ[N[x], Complex] && TrueQ[0 <= x <= 1]

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

BetaDistribution/: PDF[BetaDistribution[p_, q_],
	x_?((Negative[#] || Positive[#-1])&)] := 0 /; ParameterQ[
		BetaDistribution[p, q]]
BetaDistribution/: PDF[BetaDistribution[p_, q_], x_] :=
    x^(p-1) (1-x)^(q-1) / Beta[p, q] /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: CDF[BetaDistribution[p_, q_], x_?Negative] := 0 /;
	ParameterQ[BetaDistribution[p, q]]	
BetaDistribution/: CDF[BetaDistribution[p_, q_],
	x_?(Positive[#-1]&)] := 1 /; ParameterQ[BetaDistribution[p, q]]
BetaDistribution/: CDF[BetaDistribution[p_, q_], x_] := 
    BetaRegularized[x, p, q] /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: Mean[BetaDistribution[p_, q_]] :=
    p/(p+q) /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: Variance[BetaDistribution[p_, q_]] :=
    p q /((p+q)^2 (p+q+1)) /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: StandardDeviation[BetaDistribution[p_, q_]] :=
    Sqrt[p] Sqrt[q] /
	((p + q) Sqrt[1 + p + q]) /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: Skewness[BetaDistribution[p_, q_]] :=
    2 (q-p) Sqrt[p+q+1] /
	(Sqrt[p] Sqrt[q] (p+q+2)) /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: Kurtosis[BetaDistribution[p_, q_]] :=
    3 (p+q+1) (2 (p+q)^2 + p q (p+q-6))/
        (p q (p+q+2) (p+q+3)) /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: KurtosisExcess[BetaDistribution[p_, q_]] :=
    3 (p+q+1) (2 (p+q)^2 + p q (p+q-6))/
        (p q (p+q+2) (p+q+3)) - 3 /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/:
        CharacteristicFunction[BetaDistribution[p_, q_], t_] :=
    Hypergeometric1F1[p, p+q, I t] /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: Quantile[BetaDistribution[p_, q_], fraction_] :=
    InverseBetaRegularized[0, fraction, p, q] /; QuantileQ[fraction] &&
					ParameterQ[BetaDistribution[p, q]]						
BetaDistribution/: ExpectedValue[f_Function, BetaDistribution[p_, q_],
	opts___?OptionQ] :=
   Module[{x, integral,
	   assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[ f[x] PDF[BetaDistribution[p, q], x],
		{x, 0, 1}, Assumptions -> Join[{p > 0, q > 0}, assmp]];
	  FreeQ[integral, Integrate])
   ] /; ParameterQ[BetaDistribution[p, q]]

BetaDistribution/: ExpectedValue[f_, BetaDistribution[p_, q_], x_Symbol,
	opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f PDF[BetaDistribution[p, q], x], {x, 0, 1},
		Assumptions -> Join[{p > 0, q > 0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[BetaDistribution[p, q]]

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

BetaDistribution/: Random[BetaDistribution[p_, q_]] := 
  With[{x1 = Random[GammaDistribution[p, 1]]},
    x1/(x1 + Random[GammaDistribution[q, 1]])
  ]

BetaDistribution/: RandomArray[BetaDistribution[p_, q_], dim_] :=
  With[{x1 = RandomArray[GammaDistribution[p, 1], dim]},
    x1/(x1 + RandomArray[GammaDistribution[q, 1], dim])
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

(* ========================= Cauchy Distribution ========================== *)
(* Cauchy Distribution, J & K, Vol. 1, Chap. 16 *)

CauchyDistribution/: ParameterQ[CauchyDistribution[a_:0, b_:1]] := And[
        If[FreeQ[N[a], Complex], True,
           Message[CauchyDistribution::realparm, a]; False],
        If[FreeQ[N[b], Complex], True,
           Message[CauchyDistribution::realparm, b]; False],
        If[N[b] > 0, True,
           Message[CauchyDistribution::posscale, b]; False, True]
]

CauchyDistribution::posscale =
"The scale parameter `1` is expected to be positive."

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

CauchyDistribution/: DomainQ[CauchyDistribution[a_:0, b_:1], x_] := 
	FreeQ[N[x], Complex] 

CauchyDistribution/: Domain[CauchyDistribution[a_:0, b_:1]] :=
    Interval[{-Infinity, Infinity}] /; ParameterQ[CauchyDistribution[a, b]]

CauchyDistribution/: PDF[CauchyDistribution[a_:0, b_:1], x_] :=
    1/(b Pi (1 + ((x-a)/b)^2)) /; ParameterQ[CauchyDistribution[a, b]]

CauchyDistribution/: CDF[CauchyDistribution[a_:0, b_:1], x_] :=
    ArcTan[(x-a)/b]/Pi + 1/2 /; ParameterQ[CauchyDistribution[a, b]]

CauchyDistribution/: Mean[CauchyDistribution[___]] := Indeterminate

CauchyDistribution/: Variance[CauchyDistribution[___]] := Indeterminate

CauchyDistribution/: StandardDeviation[CauchyDistribution[___]] := Indeterminate

CauchyDistribution/: Skewness[CauchyDistribution[___]] := Indeterminate
 
CauchyDistribution/: Kurtosis[CauchyDistribution[___]] := Indeterminate

CauchyDistribution/: KurtosisExcess[CauchyDistribution[___]] := Indeterminate

CauchyDistribution/:
        CharacteristicFunction[CauchyDistribution[a_:0, b_:1], t_] :=
    Exp[I a t - b t Sign[t]] /; ParameterQ[CauchyDistribution[a, b]]

CauchyDistribution/: Quantile[CauchyDistribution[a_:0, b_:1], q_] :=
    a + b Tan[Pi (q-1/2)] /; QuantileQ[q] && ParameterQ[CauchyDistribution[a, b]]

CauchyDistribution/: ExpectedValue[f_Function,
			 CauchyDistribution[a_:0, b_:1], opts___?OptionQ] :=
  Module[{x, integral,
	 assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[ f[x] PDF[CauchyDistribution[a, b], x],
			{x, -Infinity, Infinity},
			Assumptions -> Join[{b > 0, Im[a] == 0}, assmp]];
	  FreeQ[integral, Integrate])
  ] /; ParameterQ[CauchyDistribution[a, b]]

CauchyDistribution/: ExpectedValue[f_, CauchyDistribution[a_:0, b_:1],
	 x_Symbol, opts___?OptionQ] :=
  Module[{integral,
	  assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[ f PDF[CauchyDistribution[a, b], x],
			{x, -Infinity, Infinity},
			Assumptions -> Join[{b > 0, Im[a] == 0}, assmp]];
          FreeQ[integral, Integrate])
  ] /; ParameterQ[CauchyDistribution[a, b]]

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

cauchy = Compile[{{a, _Real}, {b, _Real}, {q, _Real}}, a + b Tan[Pi (q-1/2)] ]

CauchyDistribution/: Random[CauchyDistribution[a_:0, b_:1]] :=
	cauchy[a, b, Random[]]

(* NOTE: For CauchyDistribution, RandomArray provides no speedup over Random. *)
CauchyDistribution/: RandomArray[CauchyDistribution[a_:0, b_:1], dim_] := 
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[cauchy[a, b, Random[]], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ========================== Chi-square Distribution ====================== *)
(* The chi-square distribution is defined in Statistics`NormalDistribution`. *)


(* ========================== Chi Distribution =========================== *)
(* Chi Distribution, J & K, Vol. 1, Chap. 17, p. 196 *)

ChiDistribution/: ParameterQ[ChiDistribution[n_]] := And[
        If[FreeQ[N[n], Complex], True,
           Message[ChiDistribution::realparm, n]; False],
        If[N[n] > 0, True,
           Message[ChiDistribution::posdof, n]; False, True]
]

ChiDistribution::posdof =
"The degrees of freedom parameter `1` is expected to be positive."

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

ChiDistribution/: DomainQ[ChiDistribution[___], list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
ChiDistribution/: DomainQ[ChiDistribution[___], x_] :=
	FreeQ[N[x], Complex] && TrueQ[ x >= 0 ]

ChiDistribution/: Domain[ChiDistribution[___]] := Interval[{0, Infinity}]

ChiDistribution/: PDF[ChiDistribution[n_], x_?Negative] := 0 /; ParameterQ[
							ChiDistribution[n]]
ChiDistribution/: PDF[ChiDistribution[n_], x_] :=
    x^(n-1) Exp[-x^2/2] / (2^(n/2-1) Gamma[n/2]) /; ParameterQ[
							ChiDistribution[n]]

ChiDistribution/: CDF[ChiDistribution[n_], x_?Negative] := 0 /; ParameterQ[
							ChiDistribution[n]]
ChiDistribution/: CDF[ChiDistribution[n_], x_] :=
    GammaRegularized[n/2, 0, x^2/2] /; ParameterQ[ChiDistribution[n]]

ChiDistribution/: Mean[ChiDistribution[n_]] :=
    Sqrt[2] Gamma[(n+1)/2]/Gamma[n/2] /; ParameterQ[ChiDistribution[n]]

ChiDistribution/: Variance[ChiDistribution[n_]] :=
    2 (Gamma[(n+2)/2]/Gamma[n/2] -
	         (Gamma[(n+1)/2]/Gamma[n/2])^2) /; ParameterQ[ChiDistribution[n]]

ChiDistribution/: StandardDeviation[ChiDistribution[n_]] :=
    Sqrt[2 (Gamma[(n+2)/2]/Gamma[n/2] -
		(Gamma[(n+1)/2]/Gamma[n/2])^2)] /; ParameterQ[ChiDistribution[n]]

ChiDistribution/: Skewness[ChiDistribution[n_]] :=
    (2*Gamma[(1 + n)/2]^3 - 3*Gamma[n/2]*Gamma[(1 + n)/2]*Gamma[(2 + n)/2] +
       	Gamma[n/2]^2*Gamma[(3 + n)/2])/
     	(-Gamma[(1 + n)/2]^2 + Gamma[n/2]*Gamma[(2 + n)/2])^(3/2) /; ParameterQ[
							ChiDistribution[n]]

ChiDistribution/: Kurtosis[ChiDistribution[n_]] :=
    (-3*Gamma[(1 + n)/2]^4 +
             6*Gamma[n/2]*Gamma[(1 + n)/2]^2*Gamma[(2 + n)/2] -
             4*Gamma[n/2]^2*Gamma[(1 + n)/2]*Gamma[(3 + n)/2] +
             Gamma[n/2]^3*Gamma[(4 + n)/2])/
           (Gamma[(1 + n)/2]^2 - Gamma[n/2]*Gamma[(2 + n)/2])^2 /; ParameterQ[
							ChiDistribution[n]]

ChiDistribution/: KurtosisExcess[ChiDistribution[n_]] :=
    Kurtosis[ChiDistribution[n]] - 3 /; ParameterQ[ChiDistribution[n]]

ChiDistribution/: CharacteristicFunction[ChiDistribution[n_], t_] :=
   (
   (I*2^(1/2)*t*Gamma[1/2 + n/2]*Hypergeometric1F1[1/2 + n/2, 3/2, -t^2/2])/
    Gamma[n/2] + Hypergeometric1F1[n/2, 1/2, -t^2/2]
   ) /; ParameterQ[ChiDistribution[n]]

	(*   Integrate[PDF[ChiDistribution[n], x] Exp[I t x], {x, 0, Infinity},
		Assumptions -> {Im[t] == 0}] *)

ChiDistribution/: Quantile[ChiDistribution[n_], q_] :=
    Sqrt[2] Sqrt[InverseGammaRegularized[n/2, 0, q]] /; QuantileQ[q] &&
				 ParameterQ[ChiDistribution[n]]

ChiDistribution/: ExpectedValue[f_Function, ChiDistribution[n_],
	opts___?OptionQ] :=
    Module[{x, integral,
	    assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
     (
	integral
     ) /; (integral = Integrate[ f[x] PDF[ChiDistribution[n], x],
			{x, 0, Infinity}, Assumptions -> Join[{n > 0}, assmp]];
	   FreeQ[integral, Integrate])
    ] /; ParameterQ[ChiDistribution[n]]

ChiDistribution/: ExpectedValue[f_, ChiDistribution[n_], x_Symbol,
	opts___?OptionQ] :=
   Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[ f PDF[ChiDistribution[n], x], {x, 0, Infinity},
		Assumptions -> Join[{n > 0}, assmp]];
	  FreeQ[integral, Integrate])
   ] /; ParameterQ[ChiDistribution[n]] 

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

ChiDistribution/: Random[ChiDistribution[n_]] :=
	Sqrt[Random[ChiSquareDistribution[n]]]

ChiDistribution/: RandomArray[ChiDistribution[n_], dim_] :=
   (
	Sqrt[RandomArray[ChiSquareDistribution[n], dim]]
   ) /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* =================== Noncentral Chi-Square Distribution ================ *)
(* Noncentral Chi-Square Distribution, K & J, Vol. 2, Chap. 28 *)

NoncentralChiSquareDistribution[n_, lambda_] := ChiSquareDistribution[n] /;
	lambda == 0

NoncentralChiSquareDistribution/: ParameterQ[
	NoncentralChiSquareDistribution[n_, lambda_]] := And[
   If[FreeQ[N[n], Complex], True,
           Message[NoncentralChiSquareDistribution::realparm, n]; False],
   If[FreeQ[N[lambda], Complex], True,
           Message[NoncentralChiSquareDistribution::realparm, lambda]; False],
   If[N[n] > 0, True,
      Message[NoncentralChiSquareDistribution::posdof, n]; False, True],
   If[N[lambda] > 0, True,
      Message[NoncentralChiSquareDistribution::posnoncent, lambda]; False, True]
]

NoncentralChiSquareDistribution::posdof =
"The degrees of freedom parameter `1` is expected to be positive."

NoncentralChiSquareDistribution::posnoncent =
"The noncentrality parameter `1` is expected to be positive."

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

NoncentralChiSquareDistribution/: DomainQ[
	NoncentralChiSquareDistribution[n_, lambda_], list_?VectorQ] :=
		FreeQ[N[list], Complex] &&
		Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
NoncentralChiSquareDistribution/: DomainQ[
	NoncentralChiSquareDistribution[n_, lambda_], x_] :=
		FreeQ[N[x], Complex] && TrueQ[x >= 0]

NoncentralChiSquareDistribution/:
        Domain[NoncentralChiSquareDistribution[___]] := Interval[{0, Infinity}]

NoncentralChiSquareDistribution/:
        PDF[NoncentralChiSquareDistribution[n_, lambda_], x_?Negative] := 0 /;
			ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
NoncentralChiSquareDistribution/:
        PDF[NoncentralChiSquareDistribution[n_, lambda_], x_] :=
         Exp[-(x+lambda)/2] x^(n/2-1) /(2^(n/2)) *
             Hypergeometric0F1Regularized[n/2, (lambda x)/4] /;
			ParameterQ[NoncentralChiSquareDistribution[n, lambda]]	

NoncentralChiSquareDistribution/:
        CDF[NoncentralChiSquareDistribution[n_, lambda_], x_?Negative] := 0 /;
			ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
NoncentralChiSquareDistribution/:
        CDF[NoncentralChiSquareDistribution[n_, lambda_], x_] :=
    Module[{t, nint},
    (* density = PDF[NoncentralChiSquareDistribution[n, lambda], t]; *)
        If[NumberQ[n] && NumberQ[lambda] && NumberQ[x],
	    nint = Re[
  		NIntegrate[Exp[-(t+lambda)/2] t^(n/2-1) /(2^(n/2)) *
		Hypergeometric0F1Regularized[n/2, (lambda t)/4], {t, 0, x}]];
	    If[nint > 1.0, 1.0, nint, nint],
	    t = If[SameQ[Head[x],Symbol], Unique[ToString[x]], Unique[]];
            Integrate[Exp[-(t+lambda)/2] t^(n/2-1) /(2^(n/2)) *
                Hypergeometric0F1Regularized[n/2, (lambda t)/4], {t, 0, x}] 
        ]
    ] /; ParameterQ[NoncentralChiSquareDistribution[n, lambda]]


NoncentralChiSquareDistribution/:
     Mean[NoncentralChiSquareDistribution[n_, lambda_]] :=
	n + lambda /; ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
 
NoncentralChiSquareDistribution/:
     Variance[NoncentralChiSquareDistribution[n_, lambda_]] :=
	2 n + 4 lambda /; ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
 
NoncentralChiSquareDistribution/:
 StandardDeviation[NoncentralChiSquareDistribution[n_, lambda_]] :=
  Sqrt[2 n + 4 lambda] /; ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
 
NoncentralChiSquareDistribution/:
        Skewness[NoncentralChiSquareDistribution[n_, lambda_]] :=
		Sqrt[8] (n + 3 lambda)/(n + 2 lambda)^(3/2) /;
			ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
 
NoncentralChiSquareDistribution/:
        Kurtosis[NoncentralChiSquareDistribution[n_, lambda_]] :=
		3 + 12 (n + 4 lambda)/(n + 2 lambda)^2 /;
			ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
 
NoncentralChiSquareDistribution/:
        KurtosisExcess[NoncentralChiSquareDistribution[n_, lambda_]] :=
		12 (n + 4 lambda)/(n + 2 lambda)^2 /;
			ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
 
NoncentralChiSquareDistribution/: CharacteristicFunction[
        NoncentralChiSquareDistribution[n_, lambda_], t_] :=
   Integrate[PDF[NoncentralChiSquareDistribution[n, lambda], x] Exp[I t x],
	{x, 0, Infinity}, Assumptions -> {Im[t]==0}]	/; False
 
NoncentralChiSquareDistribution/:
  Quantile[NoncentralChiSquareDistribution[n_, lambda_], q_, opts___?OptionQ] :=
	0 /; q==0 

NoncentralChiSquareDistribution/:
  Quantile[NoncentralChiSquareDistribution[n_, lambda_], q_, opts___?OptionQ] :=
        Infinity /; q==1 
 
NoncentralChiSquareDistribution/:
  Quantile[NoncentralChiSquareDistribution[n_, lambda_], q_, opts___?OptionQ] :=
    Module[{workprec, accgoal, precgoal, maxits, pdf, x,
            quantile0 = n+lambda (* initial val should be median, not mean *)	
	   },
      {workprec, accgoal, precgoal, maxits} =
	{WorkingPrecision, AccuracyGoal, PrecisionGoal, MaxIterations} /.
		{opts} /. Options[Quantile];
      {workprec, accgoal, precgoal, maxits} =
	quantileOptionCheck[workprec, accgoal, precgoal, maxits];
      pdf = Exp[-(x+lambda)/2] x^(n/2-1) /(2^(n/2)) *
			Hypergeometric0F1Regularized[n/2, (lambda x)/4];
      findQuantile[pdf, {x, 0, Null}, {0, Infinity},
	 q, quantile0, accgoal, precgoal, workprec, maxits]

    ] /; QuantileQ[q] && NumberQ[N[n]] && NumberQ[N[lambda]] &&
	 ParameterQ[NoncentralChiSquareDistribution[n, lambda]]
 
NoncentralChiSquareDistribution/: 
 ExpectedValue[f_Function, NoncentralChiSquareDistribution[n_, lambda_],
	opts___?OptionQ] :=
   Module[{x, integral,
		assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[ f[x] *
		 PDF[NoncentralChiSquareDistribution[n, lambda], x],
		 {x, 0, Infinity}, Assumptions -> Join[{n > 0, lambda > 0},
			assmp]];
	  FreeQ[integral, Integrate])
   ] /;  ParameterQ[NoncentralChiSquareDistribution[n, lambda]]

NoncentralChiSquareDistribution/:
 ExpectedValue[f_, NoncentralChiSquareDistribution[n_, lambda_], x_Symbol,
	opts___?OptionQ] :=
   Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[f *
		PDF[NoncentralChiSquareDistribution[n, lambda], x],
		{x, 0, Infinity}, Assumptions -> Join[{n > 0, lambda > 0},
			assmp]];
          FreeQ[integral, Integrate])
   ] /; ParameterQ[NoncentralChiSquareDistribution[n, lambda]]

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

NoncentralChiSquareDistribution/:
        Random[NoncentralChiSquareDistribution[n_, lambda_]] :=
   With[{delta = Sqrt[lambda/n]},
	Apply[Plus, (RandomArray[NormalDistribution[0, 1], n] + delta)^2 ]
   ] /; IntegerQ[n] && n > 0
 
NoncentralChiSquareDistribution/: (* non-integral degrees of freedom *)
        Random[NoncentralChiSquareDistribution[n_, lambda_]] :=
	Quantile[NoncentralChiSquareDistribution[n, lambda], Random[]]
 
NoncentralChiSquareDistribution/:
	RandomArray[NoncentralChiSquareDistribution[n_, lambda_], dim_] :=
  Module[{m, delta = Sqrt[lambda/n], array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = RandomArray[NormalDistribution[0, 1], {m, n}];
    array = Map[Apply[Plus, (# + delta)^2]&, array]; 
    (* Now array is a list of m NoncentralChiSquareDistribution[n, lambda]
 	variates. *)
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; IntegerQ[n] && n > 0 &&
	 (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

NoncentralChiSquareDistribution/: (* non-integral degrees of freedom *)
        RandomArray[NoncentralChiSquareDistribution[n_, lambda_], dim_] :=
  Module[{m, array},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[Random[NoncentralChiSquareDistribution[n, lambda]], {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)&]


(* ========================= Exponential Distribution ===================== *)
(* Exponential Distribution, J & K, Vol. 1, Chap. 18 *)

ExponentialDistribution/: ParameterQ[ExponentialDistribution[lambda_]] := And[
        If[FreeQ[N[lambda], Complex], True,
           Message[ExponentialDistribution::realparm, lambda]; False],
        If[N[lambda] > 0, True,
           Message[ExponentialDistribution::posscale, lambda]; False, True]
]

ExponentialDistribution::posscale =
"The parameter `1` inversely related to the scale is expected to be positive."

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

ExponentialDistribution/: DomainQ[ExponentialDistribution[lambda_],
		 list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
ExponentialDistribution/: DomainQ[ExponentialDistribution[lambda_], x_] :=
	FreeQ[N[x], Complex] && TrueQ[x >= 0]

ExponentialDistribution/: Domain[ExponentialDistribution[___]] :=
        Interval[{0, Infinity}]

ExponentialDistribution/:
	PDF[ExponentialDistribution[lambda_:1], x_?Negative] := 0 /;
				ParameterQ[ExponentialDistribution[lambda]]
ExponentialDistribution/: PDF[ExponentialDistribution[lambda_:1], x_] :=
        lambda Exp[-lambda x] /; ParameterQ[ExponentialDistribution[lambda]]

ExponentialDistribution/:
	CDF[ExponentialDistribution[lambda_:1], x_?Negative] := 0 /;
				ParameterQ[ExponentialDistribution[lambda]]
ExponentialDistribution/: CDF[ExponentialDistribution[lambda_:1], x_] :=
        1 - Exp[-lambda x] /; ParameterQ[ExponentialDistribution[lambda]]

ExponentialDistribution/: Mean[ExponentialDistribution[lambda_:1]] :=
        1/lambda /; ParameterQ[ExponentialDistribution[lambda]]

ExponentialDistribution/: Variance[ExponentialDistribution[lambda_:1]] :=
	1/lambda^2 /; ParameterQ[ExponentialDistribution[lambda]]

ExponentialDistribution/: StandardDeviation[
         ExponentialDistribution[lambda_:1]] := 1/lambda /; ParameterQ[
					ExponentialDistribution[lambda]]

ExponentialDistribution/: Skewness[ExponentialDistribution[___]] = 2

ExponentialDistribution/: Kurtosis[ExponentialDistribution[___]] = 9

ExponentialDistribution/: KurtosisExcess[ExponentialDistribution[___]] = 6

ExponentialDistribution/: CharacteristicFunction[
                              ExponentialDistribution[lambda_:1], t_] :=
        lambda/(lambda - I t) /; ParameterQ[ExponentialDistribution[lambda]]

ExponentialDistribution/: Quantile[ExponentialDistribution[lambda_:1], q_] :=
	-Log[1 - q] / lambda /; QuantileQ[q] &&  ParameterQ[
					ExponentialDistribution[lambda]]

ExponentialDistribution/: ExpectedValue[f_Function,
   		ExponentialDistribution[lambda_:1], opts___?OptionQ] :=
 Module[{x, integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
  (
	integral
  ) /; (integral = Integrate[ f[x] PDF[ExponentialDistribution[lambda], x],
			{x, 0, Infinity}, Assumptions -> Join[{lambda > 0},
			assmp]];
	FreeQ[integral, Integrate])
 ] /; ParameterQ[ExponentialDistribution[lambda]]

ExponentialDistribution/: ExpectedValue[f_, ExponentialDistribution[lambda_:1],
	x_Symbol, opts___?OptionQ] :=
 Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
  (
	integral
  ) /; (integral = Integrate[ f PDF[ExponentialDistribution[lambda], x],
		{x, 0, Infinity}, Assumptions -> Join[{lambda > 0}, assmp]];
        FreeQ[integral, Integrate])
 ] /; ParameterQ[ExponentialDistribution[lambda]]

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

exponential = Compile[{{lambda, _Real}, {q, _Real}}, -Log[q]/lambda]

ExponentialDistribution/: Random[ExponentialDistribution[lambda_:1]] :=
	exponential[lambda, Random[]]

(* NOTE: For ExponentialDistribution,
	 RandomArray provides no speedup over Random. *)
ExponentialDistribution/: RandomArray[ExponentialDistribution[lambda_], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[exponential[lambda, Random[]], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ======================= ExtremeValue Distribution ====================== *)
(* ExtremeValue (FisherTippett) Distribution, J & K, Vol. 1, Chap. 21 *)

ParameterQ[ExtremeValueDistribution[alpha_, beta_]] := And[
        If[FreeQ[N[alpha], Complex], True,
           Message[ExtremeValueDistribution::realparm, alpha]; False],
        If[FreeQ[N[beta], Complex], True,
           Message[ExtremeValueDistribution::realparm, beta]; False],
        If[N[beta] > 0, True,
           Message[ExtremeValueDistribution::posscale, beta]; False, True]
]

ExtremeValueDistribution::posscale =
"The scale parameter `1` is expected to be positive."

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

ExtremeValueDistribution/: DomainQ[ExtremeValueDistribution[alpha_, beta_],
		 x_] := FreeQ[N[x], Complex]

ExtremeValueDistribution/: Domain[ExtremeValueDistribution[___]] :=
        Interval[{-Infinity, Infinity}]

ExtremeValueDistribution/: PDF[
                ExtremeValueDistribution[alpha_, beta_], x_] :=
        Exp[-Exp[(alpha - x)/beta] + (alpha - x)/beta]/beta /;
				ParameterQ[ExtremeValueDistribution[alpha, beta]]

ExtremeValueDistribution/:
        CDF[ExtremeValueDistribution[alpha_, beta_], x_] :=
        Exp[-Exp[-(x-alpha)/beta]] /; ParameterQ[
				ExtremeValueDistribution[alpha, beta]]

ExtremeValueDistribution/: Mean[ExtremeValueDistribution[alpha_, beta_]] :=
	alpha + EulerGamma beta /; ParameterQ[
				ExtremeValueDistribution[alpha, beta]]

ExtremeValueDistribution/: Variance[
                ExtremeValueDistribution[alpha_, beta_]] :=
	Pi^2 beta^2 /6 /; ParameterQ[ExtremeValueDistribution[alpha, beta]]

ExtremeValueDistribution/: StandardDeviation[
                ExtremeValueDistribution[alpha_, beta_]] :=
	Pi beta /Sqrt[6] /; ParameterQ[ExtremeValueDistribution[alpha, beta]]

ExtremeValueDistribution/: Skewness[ExtremeValueDistribution[___]] :=
	- PolyGamma[2, 1] 6^(3/2) / Pi^3 

ExtremeValueDistribution/: Kurtosis[ExtremeValueDistribution[___]] = 27/5 

ExtremeValueDistribution/: KurtosisExcess[ExtremeValueDistribution[___]] =
	12/5 

ExtremeValueDistribution/: CharacteristicFunction[
                ExtremeValueDistribution[alpha_, beta_], t_] :=
        Exp[alpha I t] Gamma[1 - beta I t] /; ParameterQ[
					ExtremeValueDistribution[alpha, beta]]

ExtremeValueDistribution/: Quantile[
                ExtremeValueDistribution[alpha_, beta_], q_] :=
	alpha - beta Log[Log[1/q]] /; QuantileQ[q] &&
			ParameterQ[ExtremeValueDistribution[alpha, beta]]

ExtremeValueDistribution/: ExpectedValue[f_Function, 
  		ExtremeValueDistribution[alpha_, beta_], opts___?OptionQ] :=
  Module[{x, integral,
	  assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	integral = Integrate[ f[x] *
		 PDF[ExtremeValueDistribution[alpha, beta], x],
			{x, -Infinity, Infinity},
			Assumptions -> Join[{Im[alpha]==0, beta>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[ExtremeValueDistribution[alpha, beta]]

ExtremeValueDistribution/: ExpectedValue[f_,
	 ExtremeValueDistribution[alpha_, beta_], x_Symbol, opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	 integral = Integrate[ f PDF[ExtremeValueDistribution[alpha, beta], x],
		{x, -Infinity, Infinity},
		Assumptions -> Join[{Im[alpha]==0, beta>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[ExtremeValueDistribution[alpha, beta]]

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

extremevalue = Compile[{{alpha, _Real}, {beta, _Real}, {q, _Real}},
			alpha - beta Log[Log[1/q]]   ]
	
ExtremeValueDistribution/: Random[ExtremeValueDistribution[alpha_, beta_]] :=
	extremevalue[alpha, beta, Random[]]

(* NOTE: For ExtremeValueDistribution,
	 RandomArray provides no speedup over Random. *)
ExtremeValueDistribution/: RandomArray[ExtremeValueDistribution[alpha_, beta_],
	 dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[extremevalue[alpha, beta, Random[]], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ============================ F-ratio Distribution ======================== *)
(* F-ratio distribution is defined in Statistics`NormalDistribution` *)


(* ==================== Noncentral F-ratio Distribution ===================== *)
(* Noncentral F-ratio Distribution, J & K, Vol. 2, Chap. 30 *)

NoncentralFRatioDistribution[n1_, n2_, lambda_] := FRatioDistribution[n1, n2] /;
	lambda == 0

NoncentralFRatioDistribution/: ParameterQ[
	NoncentralFRatioDistribution[n1_, n2_, lambda_]] := And[
   If[FreeQ[N[n1], Complex], True,
           Message[NoncentralFRatioDistribution::realparm, n1]; False],
   If[FreeQ[N[n2], Complex], True,
           Message[NoncentralFRatioDistribution::realparm, n2]; False],
   If[FreeQ[N[lambda], Complex], True,
           Message[NoncentralFRatioDistribution::realparm, lambda]; False],
   If[N[n1] > 0, True,
      Message[NoncentralFRatioDistribution::posdof, n1]; False, True],
   If[N[n2] > 0, True,
      Message[NoncentralFRatioDistribution::posdof, n2]; False, True],
   If[N[lambda] > 0, True,
      Message[NoncentralFRatioDistribution::posnoncent, lambda]; False, True]
]

NoncentralFRatioDistribution::posdof =
"The degrees of freedom parameter `1` is expected to be positive."

NoncentralFRatioDistribution::posnoncent =
"The noncentrality parameter `1` is expected to be positive."

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

NoncentralFRatioDistribution/: DomainQ[
	NoncentralFRatioDistribution[n1_, n2_, lambda_], list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
NoncentralFRatioDistribution/: DomainQ[
	NoncentralFRatioDistribution[n1_, n2_, lambda_], x_] :=
	FreeQ[N[x], Complex] && TrueQ[x >= 0]

NoncentralFRatioDistribution/:
        Domain[NoncentralFRatioDistribution[___]] := Interval[{0, Infinity}]

NoncentralFRatioDistribution/:
     PDF[NoncentralFRatioDistribution[n1_, n2_, lambda_], x_?Negative] := 0 /;
		ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]]
NoncentralFRatioDistribution/:
     PDF[NoncentralFRatioDistribution[n1_, n2_, lambda_], x_] :=
        Exp[-lambda/2] Sqrt[x]^(n1-2) n1^(n1/2) n2^(n2/2) /
          ((n2 + x n1)^((n1+n2)/2) Beta[n1/2, n2/2]) *
          Hypergeometric1F1[(n1+n2)/2, n1/2,
          lambda n1 x/(2 (n2 + n1 x))] /; ParameterQ[
			NoncentralFRatioDistribution[n1, n2, lambda]]

NoncentralFRatioDistribution/:
     CDF[NoncentralFRatioDistribution[n1_, n2_, lambda_], x_?Negative] := 0 /;
			ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]]
NoncentralFRatioDistribution/:
        CDF[NoncentralFRatioDistribution[n1_, n2_, lambda_], x_] :=
    Module[ {t, coeff, nint},
	  (*  density = PDF[NoncentralFRatioDistribution[n1,n2,lambda], t];  *)
	  coeff = (E^(-lambda/2)*n1^(n1/2)*n2^(n2/2)*Gamma[n1/2]*
          Gamma[1 + n2/2])/(Beta[n1/2, n2/2]*Gamma[(n1 + n2)/2]);
      If[NumberQ[n1] && NumberQ[n2] && NumberQ[lambda] && NumberQ[x],  
          nint = Re[
		coeff NIntegrate[ Exp[(lambda*n1*t)/(2*(n2 + n1*t))] * 
		  t^((-2 + n1)/2)*(n2 + n1*t)^((-n1 - n2)/2) * 
		  LaguerreL[n2/2, -1 + n1/2,-(lambda*n1*t)/(2*(n2 + n1*t))],
		  {t, 0, x}]];
	  If[nint > 1.0, 1.0, nint, nint],
        (* else *)
	    t = If[SameQ[Head[x],Symbol], Unique[ToString[x]], Unique[]];
            coeff Integrate[Exp[(lambda*n1*t)/(2*(n2 + n1*t))] * 
		t^((-2 + n1)/2)*(n2 + n1*t)^((-n1 - n2)/2) * 
		  LaguerreL[n2/2, -1 + n1/2,-(lambda*n1*t)/(2*(n2 + n1*t))], 
			  {t, 0, x}]
        ]
    ] /; ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]]

NoncentralFRatioDistribution/:
        Mean[NoncentralFRatioDistribution[n1_, n2_, lambda_]] :=
	     n2(n1 + lambda)/n1/(n2 - 2)  /;
     ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]] && !TrueQ[n2 <= 2]

NoncentralFRatioDistribution/:
        Variance[NoncentralFRatioDistribution[n1_, n2_, lambda_]] :=
             2 (n2/n1)^2 ( (n1+lambda)^2 + (n1+2 lambda)(n2-2) ) /
			  (n2-2)^2 / (n2-4)  /;
    ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]] && !TrueQ[n2 <= 4]	

NoncentralFRatioDistribution/:
        StandardDeviation[NoncentralFRatioDistribution[n1_, n2_, lambda_]] :=
	    Sqrt[ Variance[NoncentralFRatioDistribution[n1, n2, lambda]] ] /;
     ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]] && !TrueQ[n2 <= 4]
 
NoncentralFRatioDistribution/:
    Skewness[NoncentralFRatioDistribution[n1_, n2_, lambda_]] :=
      Module[{mean = uncenteredmu[n1, n2, lambda, 1], var, mu3},
	     var = uncenteredmu[n1, n2, lambda, 2] - mean^2;
	     mu3 = uncenteredmu[n1, n2, lambda, 3] - 3 mean var - mean^3;
	     mu3/var^(3/2) 
      ] /; ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]] &&
								 !TrueQ[n2 <= 6]

NoncentralFRatioDistribution/:
    Kurtosis[NoncentralFRatioDistribution[n1_, n2_, lambda_]] :=
      Module[{mean = uncenteredmu[n1, n2, lambda, 1], var, mu3, mu4},
	     var = uncenteredmu[NoncentralFRatioDistribution[n1, n2, lambda],
			 2] - mean^2;
	     mu3 = uncenteredmu[NoncentralFRatioDistribution[n1, n2, lambda],
			 3] - 3 mean var - mean^3;
	     mu4 = uncenteredmu[NoncentralFRatioDistribution[n1, n2, lambda],
			 4] - 4 mean mu3 - 6 mean^2 var - mean^4;
	     mu4/var^2 
      ] /; ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]] &&
								 !TrueQ[n2 <= 8]
 
NoncentralFRatioDistribution/:
    KurtosisExcess[NoncentralFRatioDistribution[n1_, n2_, lambda_]] :=
      Module[{mean = uncenteredmu[n1, n2, lambda, 1], var, mu3, mu4},
	     var = uncenteredmu[NoncentralFRatioDistribution[n1, n2, lambda],
			 2] - mean^2;
	     mu3 = uncenteredmu[NoncentralFRatioDistribution[n1, n2, lambda],
			 3] - 3 mean var - mean^3;
	     mu4 = uncenteredmu[NoncentralFRatioDistribution[n1, n2, lambda],
			 4] - 4 mean mu3 - 6 mean^2 var - mean^4;
	     mu4/var^2 - 3 
      ] /; ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]] &&
								 !TrueQ[n2 <= 8]
 
uncenteredmu[NoncentralFRatioDistribution[nu1, nu2, lambda], r_] :=
	(nu2/nu1)^r Gamma[nu1/2+r]/Gamma[nu1/2] Gamma[nu2/2-r]/Gamma[nu2/2] *
	Sum[r!/((r-j)!j!) Gamma[nu1/2]/Gamma[nu1/2+j] (lambda/2)^j, {j, 0, r}]

NoncentralFRatioDistribution/: CharacteristicFunction[
        NoncentralFRatioDistribution[n1_, n2_, lambda_], t_] :=
  Integrate[PDF[NoncentralFRatioDistribution[n1, n2, lambda], x] Exp[I t x],
	{x, 0, Infinity}, Assumptions -> {Im[t]==0}]	/; False
 
NoncentralFRatioDistribution/:
        Quantile[NoncentralFRatioDistribution[n1_, n2_, lambda_], q_,
                opts___?OptionQ] :=
	0 /; q==0

NoncentralFRatioDistribution/:
        Quantile[NoncentralFRatioDistribution[n1_, n2_, lambda_], q_,
                opts___?OptionQ] :=
        Infinity /; q==1

NoncentralFRatioDistribution/:
        Quantile[NoncentralFRatioDistribution[n1_, n2_, lambda_], q_,
		opts___?OptionQ] :=
    Module[{workprec, accgoal, precgoal, maxits, pdf, x,
            quantile0 = n2(n1 + lambda)/n1/(n2 - 2)
		 (* initial val should be median, not mean *)
           },
      {workprec, accgoal, precgoal, maxits} =
        {WorkingPrecision, AccuracyGoal, PrecisionGoal, MaxIterations} /.
                {opts} /. Options[Quantile];
      {workprec, accgoal, precgoal, maxits} =
        quantileOptionCheck[workprec, accgoal, precgoal, maxits];
      pdf = Exp[-lambda/2] Sqrt[x]^(n1-2) n1^(n1/2) n2^(n2/2) /
          ((n2 + x n1)^((n1+n2)/2) Beta[n1/2, n2/2]) *
          Hypergeometric1F1[(n1+n2)/2, n1/2, lambda n1 x/(2 (n2 + n1 x))];
      findQuantile[pdf, {x, 0, Null}, {0, Infinity},
	 q, quantile0, accgoal, precgoal, workprec, maxits]

    ] /; QuantileQ[q] && NumberQ[N[n1]] && NumberQ[N[n2]] && n2 > 2 &&
	 NumberQ[N[lambda]] &&
         ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]]

NoncentralFRatioDistribution/: ExpectedValue[f_Function, 
	NoncentralFRatioDistribution[n1_, n2_, lambda_], opts___?OptionQ] :=
  Module[{x, integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f[x] *
	 	   PDF[NoncentralFRatioDistribution[n1, n2, lambda], x],
	 	   {x, 0, Infinity},
		   Assumptions -> Join[{n1>0, n2>0, lambda>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]]

NoncentralFRatioDistribution/: ExpectedValue[f_,
	NoncentralFRatioDistribution[n1_, n2_, lambda_], x_Symbol,
	opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f *
		 PDF[NoncentralFRatioDistribution[n1, n2, lambda], x],
		 {x, 0, Infinity},
		 Assumptions -> Join[{n1>0, n2>0, lambda>0}, assmp]];
         FreeQ[integral, Integrate])
  ] /; ParameterQ[NoncentralFRatioDistribution[n1, n2, lambda]]
	
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE: do not condition Random on ParameterQ so that Random will always
        evaluate and be efficient. *)

(* I am assuming that this is always more efficient than
	Quantile[NoncentralFRatioDistribution[n1,n2,lambda], Random[]]
*)
NoncentralFRatioDistribution/:
        Random[NoncentralFRatioDistribution[n1_, n2_, lambda_]] :=
   n2/n1 Random[NoncentralChiSquareDistribution[n1, lambda]]/
	 Random[ChiSquareDistribution[n2]]

NoncentralFRatioDistribution/:
        RandomArray[NoncentralFRatioDistribution[n1_, n2_, lambda_], dim_] :=
   n2/n1 RandomArray[NoncentralChiSquareDistribution[n1, lambda], dim]/
	 RandomArray[ChiSquareDistribution[n2], dim]


(* =========================== Gamma Distribution ======================== *)
(* Gamma Distribution, J & K, Vol. 1, Chap. 17 *)

GammaDistribution[alpha_, beta_] := ExponentialDistribution[1/beta] /; alpha==1

GammaDistribution[alpha_, beta_] := ChiSquareDistribution[2 alpha] /; beta==2

ParameterQ[GammaDistribution[alpha_, beta_]] := And[
        If[FreeQ[N[alpha], Complex], True,
           Message[GammaDistribution::realparm, alpha]; False],
        If[FreeQ[N[beta], Complex], True,
           Message[GammaDistribution::realparm, beta]; False],
        If[N[alpha] > 0, True,
           Message[GammaDistribution::posparm, alpha]; False, True],
        If[N[beta] > 0, True,
           Message[GammaDistribution::posparm, beta]; False, True]
]

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

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

GammaDistribution/: DomainQ[GammaDistribution[alpha_, beta_], list_?VectorQ] :=
		FreeQ[N[list], Complex] &&
		Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
GammaDistribution/: DomainQ[GammaDistribution[alpha_, beta_], x_] :=
		FreeQ[N[x], Complex] && TrueQ[x >= 0]

GammaDistribution/: Domain[GammaDistribution[___]] := Interval[{0, Infinity}]

GammaDistribution/: PDF[GammaDistribution[alpha_, beta_], x_?Negative] := 0 /;
			ParameterQ[GammaDistribution[alpha, beta]]
GammaDistribution/: PDF[GammaDistribution[alpha_, beta_], x_] :=
        x^(alpha-1) Exp[-x/beta] / (beta^alpha Gamma[alpha]) /;
			ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: CDF[GammaDistribution[alpha_, beta_], x_?Negative] := 0 /;
			ParameterQ[GammaDistribution[alpha, beta]]	
GammaDistribution/: CDF[GammaDistribution[alpha_, beta_], x_] :=
        GammaRegularized[alpha, 0, x/beta] /; ParameterQ[
					GammaDistribution[alpha, beta]]

GammaDistribution/: Mean[GammaDistribution[alpha_, beta_]] :=
	alpha beta /; ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: Variance[GammaDistribution[alpha_, beta_]] :=
	alpha beta^2 /; ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: StandardDeviation[GammaDistribution[alpha_, beta_]] :=
	beta Sqrt[alpha] /; ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: Skewness[GammaDistribution[alpha_, beta_]] :=
	2/Sqrt[alpha] /; ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: Kurtosis[GammaDistribution[alpha_, beta_]] :=
	6/alpha + 3 /; ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: KurtosisExcess[
                        GammaDistribution[alpha_, beta_]] :=
	6/alpha /; ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: CharacteristicFunction[
                        GammaDistribution[alpha_, beta_], t_] :=
        (1 - I beta t)^(-alpha) /; ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: Quantile[GammaDistribution[alpha_, beta_], q_] :=
	beta InverseGammaRegularized[alpha, 0, q] /;
		QuantileQ[q] && ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: ExpectedValue[f_Function,
			 GammaDistribution[alpha_, beta_], opts___?OptionQ] :=
   Module[{x, integral,
	   assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[ f[x] PDF[GammaDistribution[alpha, beta], x],
			{x, 0, Infinity},
			Assumptions -> Join[{alpha>0, beta>0}, assmp]];
	  FreeQ[integral, Integrate])
   ] /; ParameterQ[GammaDistribution[alpha, beta]]

GammaDistribution/: ExpectedValue[f_,
             GammaDistribution[alpha_, beta_], x_Symbol, opts___?OptionQ] :=
   Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
     (
	integral
     ) /; (integral = Integrate[ f PDF[GammaDistribution[alpha, beta], x],
	     {x, 0, Infinity}, Assumptions -> Join[{alpha>0, beta>0}, assmp]];
           FreeQ[integral, Integrate])
   ] /; ParameterQ[GammaDistribution[alpha, beta]]
				    
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE: do not condition Random on ParameterQ so that Random will always
        evaluate and be efficient. *)

GammaDistribution/: Random[GammaDistribution[alpha_, beta_]] :=
  Which[alpha < 1,
	  beta gsGamma[alpha, alpha/E],
	alpha < 2.5,
	  a = alpha-1;
	  beta a gkm1Gamma[(alpha - 1/(6 alpha))/a, 2/a],  
	True,
	  a = alpha-1;
	  beta a gkm2Gamma[1/Sqrt[alpha], (alpha - 1/(6 alpha))/a, 2/a] 
  ]

(* NOTE: the only savings garnered from doing
	 RandomArray[GammaDistribution[alpha, beta], n]
  rather than
	 Table[Random[GammaDistribution[alpha, beta]], {n}]
  is due to the fact that the RandomArray function computes constants
  only once. *)

GammaDistribution/: RandomArray[GammaDistribution[alpha_, beta_], dim_] :=
  Module[{n, array, r, k, a, b, c},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Which[alpha < 1,
	              r = alpha/E;
 	              beta Table[gsGamma[alpha, r], {n}],
		  alpha < 2.5,
	              a = alpha-1;  b = (alpha - 1/(6 alpha))/a;  c = 2/a;
	              beta a Table[gkm1Gamma[b, c], {n}],
		  True,
		      k = 1/Sqrt[alpha];
		      a = alpha-1;  b = (alpha - 1/(6 alpha))/a;  c = 2/a;
		      beta a Table[gkm2Gamma[k, b, c], {n}]		]; 
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /;	(IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* acceptance-rejection method for Gamma variates;
   GS algorithm of Ahrens and Dieter (1974),
	pp. 24-25, 39, "Multivariate Statistical Simulation", M.E.Johnson *)

(* NOTE: The constant "r" is passed to gsGamma in an attempt to
	 minimize overhead. *)
(* NOTE: The initialization x = 1.0, t = 1.0 is done just so that the
	V4.0 Compile is happy. *)
gsGamma = Compile[{{alpha, _Real}, {r, _Real}},
	Module[{x = 1.0, t = 1.0, q = (1 + r) Random[]},
          If[q < 1,
             x = q^(1/alpha);  t = Exp[-x],
             x = 1 - Log[1 + (1-q)/r];  t = x^(alpha-1)  ];
	  While[Random[] > t,
		q = (1 + r) Random[];
                If[q < 1,
                   x = q^(1/alpha);  t = Exp[-x],
                   x = 1 - Log[1 + (1-q)/r];  t = x^(alpha-1)  ]	];
	  x
	]	]

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* ratio of uniforms technique for Gamma variates;
   GBH algorithm of Cheng and Feast (1979),
  pp. 28-29, 40-41, "Multivariate Statistical Simulation", M.E.Johnson;
  "Some Simple Gamma Variate Generators", R.C.H.Cheng & G.M.Feast,
   Applied Statistics (1979), 28, No. 3, pp. 290-295. *)

(* NOTE: The constants "b" and "c" are passed to gkm1Gamma in an attempt to
	minimize overhead. *)
(* NOTE: The initialization w = 1.0 is done just so that the
	V4.0 Compile is happy. *)
gkm1Gamma = Compile[{{b, _Real}, {c, _Real}},
   Module[{u1 = Random[], u2 = Random[], w = 1.0},
     w = b u1/u2;
     While[!(c u2 - c - 2 + w + 1/w <= 0) && !(c Log[u2] - Log[w] + w - 1 < 0),
	   u1 = Random[];  u2 = Random[];  w = b u1/u2		];
     w
   ]	] 
	
(* NOTE: The constants "k", "b", and "c" are passed to gkm2Gamma in an attempt 
	to minimize overhead. *)
(* NOTE: The initialization u2 = 1.0, w = 1.0 is done just so that the
        V4.0 Compile is happy. *)
gkm2Gamma = Compile[{{k, _Real}, {b, _Real}, {c, _Real}},
   Module[{u1 = Random[], u2 = 1.0, w = 1.0},
     u2 = u1 + k (1 - 1.86 Random[]);
     While[!(0 < u2 < 1),
	   u1 = Random[];  u2 = u1 + k (1 - 1.86 Random[])  ];
     w = b u1/u2;
     While[!(c u2 - c - 2 + w + 1/w <= 0) &&
           !(c Log[u2] - Log[w] + w - 1 < 0),
        u1 = Random[];  u2 = u1 + k (1 - 1.86 Random[]);
	While[!(0 < u2 < 1),
	      u1 = Random[];  u2 = u1 + k (1 - 1.86 Random[])  ];
	w = b u1/u2		];
     w
   ]	]
  

(* ========================== Half-Normal Distribution ===================== *)
(* Half-Normal Distribution, J & K, Vol. 1, Chap. 13, p. 93 *)

HalfNormalDistribution/: ParameterQ[HalfNormalDistribution[theta_]] := And[
        If[FreeQ[N[theta], Complex], True,
           Message[HalfNormalDistribution::realparm, theta]; False],
        If[N[theta] > 0, True,
           Message[HalfNormalDistribution::posparm, theta]; False, True]
]

HalfNormal::posparm =
"The parameter `1` inversely related to the scale is expected to be positive."

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

HalfNormalDistribution/: DomainQ[HalfNormalDistribution[theta_],
		 list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
HalfNormalDistribution/: DomainQ[HalfNormalDistribution[theta_], x_] :=
	FreeQ[N[x], Complex] && TrueQ[x >= 0]

HalfNormalDistribution/: Domain[HalfNormalDistribution[___]] :=
	 Interval[{0, Infinity}]

HalfNormalDistribution/: PDF[HalfNormalDistribution[theta_], x_?Negative] :=
				0 /; ParameterQ[HalfNormalDistribution[theta]]
HalfNormalDistribution/: PDF[HalfNormalDistribution[theta_], x_] :=
        2 theta Exp[- theta^2 x^2 / Pi] / Pi /; ParameterQ[
						HalfNormalDistribution[theta]]

HalfNormalDistribution/: CDF[HalfNormalDistribution[theta_], x_?Negative] := 0
HalfNormalDistribution/: CDF[HalfNormalDistribution[theta_], x_] :=
        Erf[theta x / Sqrt[Pi]] /; ParameterQ[HalfNormalDistribution[theta]]

HalfNormalDistribution/: Mean[HalfNormalDistribution[theta_]] :=
	1/theta /; ParameterQ[HalfNormalDistribution[theta]]

HalfNormalDistribution/: Variance[
                HalfNormalDistribution[theta_]] :=
	(Pi-2) / (2 theta^2) /; ParameterQ[HalfNormalDistribution[theta]]

HalfNormalDistribution/: StandardDeviation[
                HalfNormalDistribution[theta_]] :=
	Sqrt[Pi/2-1]/theta /; ParameterQ[HalfNormalDistribution[theta]]

HalfNormalDistribution/: Skewness[HalfNormalDistribution[___]] :=
        Sqrt[2] (4 - Pi) / Sqrt[Pi - 2]^3

HalfNormalDistribution/: Kurtosis[HalfNormalDistribution[___]] :=
 	(3 Pi^2 - 4 Pi - 12) / (Pi - 2)^2

HalfNormalDistribution/: KurtosisExcess[HalfNormalDistribution[___]] :=
        (8 Pi - 24)/(Pi - 2)^2

HalfNormalDistribution/: CharacteristicFunction[
				HalfNormalDistribution[theta_], t_] :=
   (
	(1 + I*Erfi[(Pi^(1/2)*t)/(2*theta)])/E^((Pi*t^2)/(4*theta^2))
   ) /; ParameterQ[HalfNormalDistribution[theta]]
  
   (* Integrate[PDF[HalfNormalDistribution[theta], x] Exp[I x t],
	{x, 0, Infinity}, Assumptions -> {Im[t]==0, theta > 0}] *)

HalfNormalDistribution/: Quantile[HalfNormalDistribution[theta_], q_] :=
	Sqrt[Pi] InverseErf[q] / theta /; QuantileQ[q] && ParameterQ[
					HalfNormalDistribution[theta]]

HalfNormalDistribution/: ExpectedValue[f_Function,
		HalfNormalDistribution[theta_], opts___?OptionQ] :=
  Module[{x, integral,
	  assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f[x] PDF[HalfNormalDistribution[theta], x],
			{x, 0, Infinity},
			Assumptions -> Join[{theta>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[HalfNormalDistribution[theta]]

HalfNormalDistribution/: ExpectedValue[f_, HalfNormalDistribution[theta_],
	x_Symbol, opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f PDF[HalfNormalDistribution[theta], x],
		{x, 0, Infinity}, Assumptions -> Join[{theta>0}, assmp]];
         FreeQ[integral, Integrate])
  ] /; ParameterQ[HalfNormalDistribution[theta]]

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

halfnormal = Compile[{{theta, _Real}, {q1, _Real}, {q2, _Real}},
	Sqrt[Pi/2]/theta Sqrt[-2 Log[q1]] Abs[Cos[2Pi q2]]		]

HalfNormalDistribution/: Random[HalfNormalDistribution[theta_]] :=
	halfnormal[theta, Random[], Random[]]

halfnormalpair = Compile[{{theta, _Real}, {q1, _Real}, {q2, _Real}},
   Sqrt[Pi/2]/theta Sqrt[-2 Log[q1]] {Abs[Cos[2Pi q2]], Abs[Sin[2Pi q2]]}  ]

HalfNormalDistribution/: RandomArray[HalfNormalDistribution[theta_], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Flatten[Table[halfnormalpair[theta, Random[], Random[]],
                 {Quotient[n, 2]}]];
    If[OddQ[n],
       AppendTo[array, halfnormal[theta, Random[], Random[]] ]  ];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ================= Laplace (Double Exponential) Distribution ============== *)
(* Laplace (Double Exponential) Distribution, J & K, Vol. 2, Chap. 23 *)

LaplaceDistribution/: ParameterQ[LaplaceDistribution[mu_, beta_]] := And[
        If[FreeQ[N[mu], Complex], True,
           Message[LaplaceDistribution::realparm, mu]; False],
        If[FreeQ[N[beta], Complex], True,
           Message[LaplaceDistribution::realparm, beta]; False],
        If[N[beta] > 0, True,
           Message[LaplaceDistribution::posscale, beta]; False, True]
]

LaplaceDistribution::posscale =
"The scale parameter `1` is expected to be positive."

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

LaplaceDistribution/: DomainQ[LaplaceDistribution[mu_, beta_], x_] :=
	FreeQ[N[x], Complex]

LaplaceDistribution/: Domain[LaplaceDistribution[___]] :=
        Interval[{-Infinity, Infinity}]

LaplaceDistribution/: PDF[LaplaceDistribution[mu_, beta_], x_] :=
    (Exp[-Sign[x-mu] (x-mu)/beta] / (2 beta)) /; ParameterQ[
					LaplaceDistribution[mu, beta]]

LaplaceDistribution/: CDF[LaplaceDistribution[mu_, beta_], x_] :=
    ((1 + Sign[x-mu] (1 - Exp[-Sign[x-mu] (x-mu)/beta])) / 2) /; ParameterQ[
					LaplaceDistribution[mu, beta]]

LaplaceDistribution/: Mean[LaplaceDistribution[mu_, beta_]] :=
	mu /; ParameterQ[LaplaceDistribution[mu, beta]]

LaplaceDistribution/: Variance[LaplaceDistribution[mu_, beta_]] :=
	2 beta^2 /; ParameterQ[LaplaceDistribution[mu, beta]]

LaplaceDistribution/: StandardDeviation[LaplaceDistribution[mu_, beta_]] :=
	Sqrt[2] beta /; ParameterQ[LaplaceDistribution[mu, beta]]

LaplaceDistribution/: Skewness[LaplaceDistribution[___]] := 0

LaplaceDistribution/: Kurtosis[LaplaceDistribution[___]] := 6 

LaplaceDistribution/: KurtosisExcess[LaplaceDistribution[___]] := 3 

LaplaceDistribution/: CharacteristicFunction[
                          LaplaceDistribution[mu_, beta_], t_] :=
      Exp[I mu t]/(1 + beta^2 t^2) /; ParameterQ[LaplaceDistribution[mu, beta]]

LaplaceDistribution/: Quantile[LaplaceDistribution[mu_, beta_], q_] :=
	mu - beta Sign[2 q - 1] Log[1 - Sign[2 q - 1]
		 (2 q - 1)] /; QuantileQ[q] && ParameterQ[
					LaplaceDistribution[mu, beta]]

LaplaceDistribution/: ExpectedValue[f_Function,
	 LaplaceDistribution[mu_, beta_], opts___?OptionQ] :=
  Module[{x, integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f[x] PDF[LaplaceDistribution[mu, beta], x],
			{x, -Infinity, Infinity},
			Assumptions -> Join[{Im[mu]==0, beta>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[LaplaceDistribution[mu, beta]]

LaplaceDistribution/: ExpectedValue[f_, LaplaceDistribution[mu_, beta_],
	 x_Symbol, opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[ f PDF[LaplaceDistribution[mu, beta], x],
		{x, -Infinity, Infinity},
		Assumptions -> Join[{Im[mu]==0, beta>0}, assmp]];
          FreeQ[integral, Integrate])
  ] /; ParameterQ[LaplaceDistribution[mu, beta]]

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

laplace = Compile[{{mu, _Real}, {beta, _Real}, {q, _Real}},
		mu - beta Sign[2 q - 1] Log[1 - Sign[2 q - 1] (2 q - 1)]  ]

LaplaceDistribution/: Random[LaplaceDistribution[mu_, beta_]] :=
	laplace[mu, beta, Random[]]

(* NOTE: For LaplaceDistribution,
	 RandomArray provides no speedup over Random. *)
LaplaceDistribution/: RandomArray[LaplaceDistribution[mu_, beta_], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[laplace[mu, beta, Random[]], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ========================= Log-Normal Distribution ======================= *)
(* Log-Normal Distribution, J & K, Vol. 1, Chap. 14, p. 112 *)

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

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

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

LogNormalDistribution/: DomainQ[LogNormalDistribution[mu_, sigma_],
		 list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
LogNormalDistribution/: DomainQ[LogNormalDistribution[mu_, sigma_], x_] :=
	FreeQ[N[x], Complex] && TrueQ[x >= 0]

LogNormalDistribution/: Domain[LogNormalDistribution[___]] :=
	 Interval[{0, Infinity}]

LogNormalDistribution/:
	PDF[LogNormalDistribution[mu_, sigma_],
		 x_?((# <= 0)&)] := 0 /;
				ParameterQ[LogNormalDistribution[mu, sigma]]
LogNormalDistribution/: PDF[LogNormalDistribution[mu_, sigma_], x_] :=
      ( Exp[-(Log[x] - mu)^2 / (2 sigma^2)] /
				 (Sqrt[2 Pi] sigma x) ) /; ParameterQ[
					LogNormalDistribution[mu, sigma]]

LogNormalDistribution/:
	CDF[LogNormalDistribution[mu_, sigma_], x_?Negative] := 0 /;
				ParameterQ[LogNormalDistribution[mu, sigma]]
LogNormalDistribution/: CDF[LogNormalDistribution[mu_, sigma_], x_] :=
      (Erf[(-mu + Log[x]) / (Sqrt[2] sigma)] + 1)/2 /; ParameterQ[
					LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: Mean[LogNormalDistribution[mu_, sigma_]] :=
      Exp[mu + sigma^2/2] /; ParameterQ[LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: Variance[LogNormalDistribution[mu_, sigma_]] :=
      Exp[2 mu + sigma^2] (Exp[sigma^2] - 1) /; ParameterQ[
					LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: StandardDeviation[
                LogNormalDistribution[mu_, sigma_]] :=
      Sqrt[Exp[2 mu + sigma^2] (Exp[sigma^2] - 1)] /; ParameterQ[
					LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: Skewness[LogNormalDistribution[mu_, sigma_]] :=
      (Exp[sigma^2] + 2) Sqrt[Exp[sigma^2] - 1] /; ParameterQ[
					LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: Kurtosis[LogNormalDistribution[mu_, sigma_]] :=
      Exp[4 sigma^2] + 2 Exp[3 sigma^2] + 3 Exp[2 sigma^2] /; ParameterQ[
					LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: KurtosisExcess[LogNormalDistribution[mu_, sigma_]] :=
   Exp[4 sigma^2] + 2 Exp[3 sigma^2] + 3 Exp[2 sigma^2] - 3 /; ParameterQ[
					LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: CharacteristicFunction[
                LogNormalDistribution[mu_, sigma_], t_] :=
   Integrate[PDF[LogNormalDistribution[mu, sigma], x] Exp[I x t],
	{x, 0, Infinity}, Assumptions -> {Im[t]==0}]	/; False

LogNormalDistribution/: Quantile[LogNormalDistribution[mu_, sigma_], q_] :=
   Exp[InverseErf[2 q - 1] Sqrt[2] sigma + mu] /; QuantileQ[q] &&
				ParameterQ[LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: ExpectedValue[f_Function, 
		LogNormalDistribution[mu_, sigma_], opts___?OptionQ] :=
  Module[{x, integral, 
	  assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f[x] *
			PDF[LogNormalDistribution[mu, sigma], x],
			{x, 0, Infinity},
			Assumptions -> Join[{Im[mu]==0, sigma>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[LogNormalDistribution[mu, sigma]]

LogNormalDistribution/: ExpectedValue[f_, LogNormalDistribution[mu_, sigma_],
	 x_Symbol, opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f PDF[LogNormalDistribution[mu, sigma], x],
		{x, 0, Infinity},
		Assumptions -> Join[{Im[mu]==0, sigma>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[LogNormalDistribution[mu, sigma]]

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

lognormal = Compile[{{mu, _Real}, {sigma, _Real}, {q1, _Real}, {q2, _Real}},
		Exp[mu + sigma Sqrt[-2 Log[q1]] Cos[2Pi q2]]	]

LogNormalDistribution/: Random[LogNormalDistribution[mu_, sigma_]] :=
	lognormal[mu, sigma, Random[], Random[]]

lognormalpair = Compile[{{mu, _Real}, {sigma, _Real}, {q1, _Real}, {q2, _Real}},
   With[{k = sigma Sqrt[-2 Log[q1]]},
	{Exp[mu + k Cos[2Pi q2]], Exp[mu + k Sin[2Pi q2]]}
   ]	]

LogNormalDistribution/: RandomArray[LogNormalDistribution[mu_, sigma_], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Flatten[Table[lognormalpair[mu, sigma, Random[], Random[]],
                 {Quotient[n, 2]}]];
    If[OddQ[n],
       AppendTo[array, lognormal[mu, sigma, Random[], Random[]] ]  ];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ========================== Logistic Distribution ======================== *)
(* Logistic Distribution, K & J, Vol. 2, Chap. 22, p. 1 *)

LogisticDistribution/: ParameterQ[LogisticDistribution[mu_, beta_]] := And[
        If[FreeQ[N[mu], Complex], True,
           Message[LogisticDistribution::realparm, mu]; False],
        If[FreeQ[N[beta], Complex], True,
           Message[LogisticDistribution::realparm, beta]; False],
        If[N[beta] > 0, True,
           Message[LogisticDistribution::posscale, beta]; False, True]
]

LogisticDistribution::posscale =
"The scale parameter `1` is expected to be positive."

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

LogisticDistribution/: DomainQ[LogisticDistribution[mu_, beta_], x_] :=
	FreeQ[N[x], Complex]

LogisticDistribution/: Domain[LogisticDistribution[___]] :=
        Interval[{-Infinity, Infinity}]

LogisticDistribution/: PDF[LogisticDistribution[mu_:0, beta_:1], x_] :=
 	Exp[-(x-mu)/beta] / (beta (1 + Exp[-(x-mu)/beta])^2) /; ParameterQ[  
			LogisticDistribution[mu, beta]]

LogisticDistribution/: CDF[LogisticDistribution[mu_:0, beta_:1], x_] :=
        1/(1 + Exp[-(x - mu)/beta]) /; ParameterQ[
                        LogisticDistribution[mu, beta]]

LogisticDistribution/: Mean[LogisticDistribution[mu_:0, beta_]] := 
	mu  /; ParameterQ[LogisticDistribution[mu, beta]]

LogisticDistribution/: Variance[LogisticDistribution[mu_, beta_:1]] :=
	Pi^2 beta^2 / 3 /; 	ParameterQ[LogisticDistribution[mu, beta]]

LogisticDistribution/: StandardDeviation[LogisticDistribution[mu_, beta_:1]] :=
	Pi beta / Sqrt[3] /; ParameterQ[LogisticDistribution[mu, beta]]

LogisticDistribution/: Skewness[LogisticDistribution[___]] := 0

LogisticDistribution/: Kurtosis[LogisticDistribution[___]] := 21/5 

LogisticDistribution/: KurtosisExcess[LogisticDistribution[___]] := 6/5 

LogisticDistribution/: CharacteristicFunction[
                          LogisticDistribution[mu_:0, beta_:1], t_] :=
        I Exp[I mu t] Pi beta t / Sin[I Pi beta t] /; ParameterQ[
					LogisticDistribution[mu, beta]]

LogisticDistribution/: Quantile[LogisticDistribution[mu_:0, beta_:1], q_] :=
	mu - beta Sign[beta] Log[1-q] +
		 beta Sign[beta] Log[q] /; QuantileQ[q] && ParameterQ[
					LogisticDistribution[mu, beta]]

LogisticDistribution/: ExpectedValue[f_Function,
		LogisticDistribution[mu_:0, beta_:1], opts___?OptionQ] :=
  Module[{x, integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	 integral = Integrate[ f[x] PDF[LogisticDistribution[mu, beta], x],
		{x, -Infinity, Infinity},
	        Assumptions -> Join[{beta > 0, Im[mu] == 0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[LogisticDistribution[mu, beta]]

LogisticDistribution/: ExpectedValue[f_, LogisticDistribution[mu_:0, beta_:1],
	 x_Symbol, opts___?OptionQ] :=
 Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
  (
	integral
  ) /; (
	integral = Integrate[ f PDF[LogisticDistribution[mu, beta], x],
		{x, -Infinity, Infinity},
                Assumptions -> Join[{beta > 0, Im[mu] == 0}, assmp]];
        FreeQ[integral, Integrate])
 ] /; ParameterQ[LogisticDistribution[mu, beta]]

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

logistic = Compile[{{mu, _Real}, {beta, _Real}, {q, _Real}},
		mu - beta Sign[beta] Log[1-q] + beta Sign[beta] Log[q] ]

LogisticDistribution/: Random[LogisticDistribution[mu_:0, beta_:1]] :=
		logistic[mu, beta, Random[]]

(* NOTE: For LogisticDistribution,
	 RandomArray provides no speedup over Random. *)
LogisticDistribution/: RandomArray[LogisticDistribution[mu_:0, beta_:1],
	 dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[logistic[mu, beta, Random[]], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ========================== Normal Distribution ======================== *)
(* The normal distribution is defined in Statistics`NormalDistribution` *)


(* ========================== Rayleigh Distribution ======================= *)
(* Rayleigh Distribution, K & J, Vol. 1, Chap. 17, p. 197 *)

ParameterQ[RayleighDistribution[sigma_]] := And[
        If[FreeQ[N[sigma], Complex], True,
           Message[RayleighDistribution::realparm, sigma]; False],
        If[N[sigma] > 0, True,
           Message[RayleighDistribution::posparm, sigma]; False, True]
]

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

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

RayleighDistribution/: DomainQ[RayleighDistribution[sigma_], list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
RayleighDistribution/: DomainQ[RayleighDistribution[sigma_], x_] :=
	FreeQ[N[x], Complex] && TrueQ[x >= 0]

RayleighDistribution/: Domain[RayleighDistribution[___]] :=
        Interval[{0, Infinity}]

RayleighDistribution/: PDF[RayleighDistribution[sigma_:1], x_?Negative] := 0 /;
				ParameterQ[RayleighDistribution[sigma]]
RayleighDistribution/: PDF[RayleighDistribution[sigma_:1], x_] :=
        x Exp[-x^2/(2 sigma^2)]/sigma^2 /; ParameterQ[
					RayleighDistribution[sigma]]

RayleighDistribution/: CDF[RayleighDistribution[sigma_:1], x_?Negative] := 0 /;
				ParameterQ[RayleighDistribution[sigma]]
RayleighDistribution/: CDF[RayleighDistribution[sigma_:1], x_] :=
        1 - Exp[- x^2 / (2 sigma^2)] /; ParameterQ[RayleighDistribution[sigma]]

RayleighDistribution/: Mean[RayleighDistribution[sigma_:1]] :=
	sigma Sqrt[Pi/2] /; ParameterQ[RayleighDistribution[sigma]]

RayleighDistribution/: Variance[RayleighDistribution[sigma_:1]] :=
	sigma^2 (2 - Pi/2) /; ParameterQ[RayleighDistribution[sigma]]

RayleighDistribution/: StandardDeviation[RayleighDistribution[sigma_:1]] :=
	sigma Sqrt[2 - Pi/2] /; ParameterQ[RayleighDistribution[sigma]]

RayleighDistribution/: Skewness[RayleighDistribution[___]] :=
        (Pi - 3) Sqrt[Pi/2] / (2 - Pi/2)^(3/2)

RayleighDistribution/: Kurtosis[RayleighDistribution[___]] :=
        (32 - 3 Pi^2)/(4 - Pi)^2

RayleighDistribution/: KurtosisExcess[RayleighDistribution[___]] =
        (32 - 3 Pi^2)/(4 - Pi)^2 - 3
 
RayleighDistribution/: CharacteristicFunction[
                           RayleighDistribution[sigma_:1], t_] :=
	1 + I t sigma Sqrt[Pi/2] Exp[-t^2 sigma^2/2] *
		Erf[-I t sigma/Sqrt[2],Infinity] /; ParameterQ[
						RayleighDistribution[sigma]]
 
RayleighDistribution/: Quantile[RayleighDistribution[sigma_:1], q_] :=
	sigma Sqrt[Log[1/(1-q)^2]] /; QuantileQ[q] && ParameterQ[
					RayleighDistribution[sigma]]
 
RayleighDistribution/: ExpectedValue[f_Function,
		 RayleighDistribution[sigma_:1], opts___?OptionQ] :=
  Module[{x, integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f[x] PDF[RayleighDistribution[sigma], x],
			{x, 0, Infinity},
			Assumptions -> Join[{sigma > 0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[RayleighDistribution[sigma]]

RayleighDistribution/: ExpectedValue[f_, RayleighDistribution[sigma_:1],
	x_Symbol, opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
    (
	integral
    ) /; (integral = Integrate[ f PDF[RayleighDistribution[sigma], x],
		{x, 0, Infinity}, Assumptions -> Join[{sigma > 0}, assmp]];
          FreeQ[integral, Integrate])
  ] /; ParameterQ[RayleighDistribution[sigma]]

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

rayleigh = Compile[{{sigma, _Real}, {q, _Real}},
			sigma Sqrt[Log[1/q^2]]  ]

RayleighDistribution/: Random[RayleighDistribution[sigma_:1]] :=
	rayleigh[sigma, Random[]]	

(* NOTE: For RayleighDistribution,
	 RandomArray provides no speedup over Random. *)
RayleighDistribution/: RandomArray[RayleighDistribution[sigma_:1], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[rayleigh[sigma, Random[]], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ======================== Student's t Distribution ======================= *)
(* Student's t distribution is defined in Statistics`NormalDistribution`. *)


(* ==================== Noncentral Student's t Distribution ================= *)
(* Noncentral Student's t Distribution, K & J, Vol. 2, Chap. 31 *)

NoncentralStudentTDistribution[n_, delta_] := StudentTDistribution[n] /;
	delta == 0

ParameterQ[NoncentralStudentTDistribution[n_, delta_]] := And[
   If[FreeQ[N[n], Complex], True,
           Message[NoncentralStudentTDistribution::realparm, n]; False],
   If[FreeQ[N[delta], Complex], True,
           Message[NoncentralStudentTDistribution::realparm, delta]; False],
   If[N[n] > 0, True,
      Message[NoncentralStudentTDistribution::posdof, n]; False, True]
]

NoncentralStudentTDistribution::posdof =
"The degrees of freedom parameter `1` is expected to be positive."

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

NoncentralStudentTDistribution/:
        Domain[NoncentralStudentTDistribution[___]] := 
		Interval[{-Infinity, Infinity}]
 
NoncentralStudentTDistribution/:
        PDF[NoncentralStudentTDistribution[n_, delta_], x_] :=
    Module[{com1, com2},
        com1 = n + x^2;
        com2 = delta^2 x^2 / (2 com1);
        n! n^(n/2) / (2^n Gamma[n/2] Exp[delta^2/2]) (  
            Sqrt[2] delta x com1^(-1 - n/2) / Gamma[(1 + n)/2] *
                Hypergeometric1F1[1 + n/2, 3/2, com2] +   
            (n + x^2)^(-1/2 - n/2) / Gamma[1 + n/2] * 
                Hypergeometric1F1[(1 + n)/2, 1/2, com2] )
    ] /; ParameterQ[NoncentralStudentTDistribution[n, delta]]
 
NoncentralStudentTDistribution/:
        CDF[NoncentralStudentTDistribution[n_, delta_], x_] :=
    Module[{t, com1, com2, coeff},
	(* density = PDF[NoncentralStudentTDistribution[n, delta], t]; *)
	com1 = n + t^2;
    com2 = delta^2 t^2 / (2 com1);
	coeff = n! n^(n/2) / (2^n Gamma[n/2] Gamma[(1+n)/2] Exp[delta^2/2]);
	If[NumberQ[n] && NumberQ[delta] && NumberQ[x],
        	nint = Re[
		  coeff NIntegrate[ 2^(1/2)*com1^(-1 - n/2)*delta*t* 
		  Hypergeometric1F1[1 + n/2, 3/2, com2] + E^com2*Pi^(1/2)* 
		  (n + t^2)^(-1/2 - n/2)*LaguerreL[n/2, -1/2, -com2], 
		          {t, -Infinity, x}]];
		If[nint > 1.0, 1.0, nint, nint],
	    t = If[SameQ[Head[x],Symbol], Unique[ToString[x]], Unique[]];
            coeff Integrate[ 2^(1/2)*com1^(-1 - n/2)*delta*t* 
			Hypergeometric1F1[1 + n/2, 3/2, com2] + E^com2*Pi^(1/2)* 
			(n + t^2)^(-1/2 - n/2)*LaguerreL[n/2, -1/2, -com2], 
		          {t, -Infinity, x}]
        ]	
    ] /; ParameterQ[NoncentralStudentTDistribution[n, delta]]

NoncentralStudentTDistribution/:
        Mean[NoncentralStudentTDistribution[n_, delta_]] :=
	   Sqrt[n/2] Gamma[(n-1)/2] / Gamma[n/2] delta  /; ParameterQ[
		NoncentralStudentTDistribution[n, delta]] && !TrueQ[n <= 1]

NoncentralStudentTDistribution/:
        Variance[NoncentralStudentTDistribution[n_, delta_]] :=
	  With[{mean = Sqrt[n/2] Gamma[(n-1)/2] / Gamma[n/2] delta},
		  n/(n-2) (1 + delta^2) - mean^2  
	  ] /; ParameterQ[NoncentralStudentTDistribution[n, delta]] &&
							 !TrueQ[n <= 2]
 
NoncentralStudentTDistribution/:
        StandardDeviation[NoncentralStudentTDistribution[n_, delta_]] :=
	Sqrt[ Variance[NoncentralStudentTDistribution[n, delta]] ] /;
	ParameterQ[NoncentralStudentTDistribution[n, delta]] && !TrueQ[n <= 2]
 
NoncentralStudentTDistribution/:
        Skewness[NoncentralStudentTDistribution[n_, delta_]] :=
	  Module[{mean = Sqrt[n/2] Gamma[(n-1)/2] / Gamma[n/2] delta,
		  variance, mu3},
	    variance = n/(n-2) (1 + delta^2) - mean^2;
	    mu3 = mean (  n(2 n-3+delta^2)/(n-2)/(n-3) - 2 variance  );
	    mu3 / variance^(3/2)  
	  ] /; ParameterQ[NoncentralStudentTDistribution[n, delta]] &&
								 !TrueQ[n <= 3] 
 
NoncentralStudentTDistribution/:
        Kurtosis[NoncentralStudentTDistribution[n_, delta_]] :=
	  Module[{mean = Sqrt[n/2] Gamma[(n-1)/2] / Gamma[n/2] delta,
	          variance, mu4},
	    variance = n/(n-2) (1 + delta^2) - mean^2;
	    mu4 = n^2 (3 + 6 delta^2 + delta^4)/(n-2)/(n-4) -
		mean^2 ( n((n+1)delta^2 + 3(3 n - 5))/(n-2)/(n-3) - 3 variance);
	    mu4 / variance^2 
	  ] /; ParameterQ[NoncentralStudentTDistribution[n, delta]] &&
								 !TrueQ[n <= 4]
 
NoncentralStudentTDistribution/:
        KurtosisExcess[NoncentralStudentTDistribution[n_, delta_]] :=
	  Module[{mean = Sqrt[n/2] Gamma[(n-1)/2] / Gamma[n/2] delta,
                  variance, mu4},
            variance = n/(n-2) (1 + delta^2) - mean^2;
            mu4 = n^2 (3 + 6 delta^2 + delta^4)/(n-2)/(n-4) -
                mean^2 ( n((n+1)delta^2 + 3(3 n - 5))/(n-2)/(n-3) - 3 variance);
            mu4 / variance^2 - 3 
          ] /; ParameterQ[NoncentralStudentTDistribution[n, delta]] &&
								 !TrueQ[n <= 4]
 
NoncentralStudentTDistribution/: CharacteristicFunction[
        NoncentralStudentTDistribution[n_, delta_], t_] :=
  Integrate[PDF[NoncentralStudentTDistribution[n, delta], x] Exp[I t x],
	{x, -Infinity, Infinity}, Assumptions -> {Im[t]==0}]	/; False

NoncentralStudentTDistribution/:
        Quantile[NoncentralStudentTDistribution[n_, delta_], q_,
		 opts___?OptionQ] := -Infinity /; q==0
  
NoncentralStudentTDistribution/:
        Quantile[NoncentralStudentTDistribution[n_, delta_], q_,
		 opts___?OptionQ] := Infinity /; q==1
 
NoncentralStudentTDistribution/:
        Quantile[NoncentralStudentTDistribution[n_, delta_], q_,
		 opts___?OptionQ] :=
    Module[{workprec, accgoal, precgoal, maxits, p0, com1, com2, pdf, x, 
            median  (* see eq. (14.1) in Johnson & Kotz
		       for an approximation to the quantile *) 
           },
      {workprec, accgoal, precgoal, maxits} =
        {WorkingPrecision, AccuracyGoal, PrecisionGoal, MaxIterations} /.
                {opts} /. Options[Quantile];
      {workprec, accgoal, precgoal, maxits} =
        quantileOptionCheck[workprec, accgoal, precgoal, maxits];
      p0 = CDF[NormalDistribution[0, 1], -delta];
      com1 = n + x^2;
      com2 = delta^2 x^2 / (2 com1);
      (* NOTE: Exp[delta^2/2] can be a very large constant! *)
      pdf = n! n^(n/2) / (2^n Gamma[n/2] Exp[delta^2/2]) (
            Sqrt[2] delta x com1^(-1 - n/2) / Gamma[(1 + n)/2] *
                Hypergeometric1F1[1 + n/2, 3/2, com2] +
            (n + x^2)^(-1/2 - n/2) / Gamma[1 + n/2] *
                Hypergeometric1F1[(1 + n)/2, 1/2, com2] );
      median = delta / ( Sqrt[2/n] Gamma[(n+1)/2] / Gamma[n/2] );
      If[p0 < q, 
         findQuantile[pdf, {x, 0, Null}, {0, Infinity}, q-p0,
		 	median, accgoal, precgoal, workprec, maxits],
         findQuantile[pdf, {x, Null, 0}, {-Infinity, 0}, p0-q,
		 	median, accgoal, precgoal, workprec, maxits]
	 
      ]

    ] /; QuantileQ[q] && NumberQ[N[n]] && NumberQ[N[delta]] &&
         ParameterQ[NoncentralStudentTDistribution[n, delta]]
 
NoncentralStudentTDistribution/: ExpectedValue[f_Function, 
		NoncentralStudentTDistribution[n_, delta_], opts___?OptionQ] :=
  Module[{x, integral,
	  assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	integral = Integrate[ f[x] *
			PDF[NoncentralStudentTDistribution[n, delta], x],
			{x, -Infinity, Infinity},
			Assumptions -> Join[{n > 0, Im[delta]==0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[NoncentralStudentTDistribution[n, delta]]

NoncentralStudentTDistribution/: ExpectedValue[f_,
   NoncentralStudentTDistribution[n_, delta_], x_Symbol, opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	integral = Integrate[ f *
		 PDF[NoncentralStudentTDistribution[n, delta], x],
		 {x, -Infinity, Infinity}, 
		 Assumptions -> Join[{n > 0, Im[delta]==0}, assmp]];
         FreeQ[integral, Integrate])
  ] /; ParameterQ[NoncentralStudentTDistribution[n, delta]]

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

NoncentralStudentTDistribution/:
        Random[NoncentralStudentTDistribution[n_, delta_]] :=
  Sqrt[n] (Random[NormalDistribution[0, 1]] + delta) /
		Random[ChiSquareDistribution[n]]	
 
NoncentralStudentTDistribution/:
	RandomArray[NoncentralStudentTDistribution[n_, delta_], dim_] :=
  (
  Sqrt[n] (RandomArray[NormalDistribution[0, 1], dim] + delta) /
		RandomArray[ChiSquareDistribution[n], dim]
  ) /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* =========================== Pareto Distribution ======================= *)
(* Pareto Distribution, K & J, Vol. 1, Chap. 19 *)

ParetoDistribution/: ParameterQ[ParetoDistribution[k_, alpha_]] := And[
        If[FreeQ[N[k], Complex], True,
           Message[ParetoDistribution::realparm, k]; False],
        If[FreeQ[N[alpha], Complex], True,
           Message[ParetoDistribution::realparm, alpha]; False],
        If[N[k] > 0, True,
           Message[ParetoDistribution::posmin, k]; False, True],
        If[N[alpha] > 0, True,
           Message[ParetoDistribution::posshape, alpha]; False, True]
]

ParetoDistribution::posmin =
"The minimum value parameter `1` is expected to be positive."

ParetoDistribution::posshape =
"The shape parameter `1` is expected to be positive."

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

ParetoDistribution/: DomainQ[ParetoDistribution[k_, alpha_], list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[# >= k], Return[False]]&, list] =!= False
ParetoDistribution/: DomainQ[ParetoDistribution[k_, alpha_], x_] :=
	FreeQ[N[x], Complex] && TrueQ[x >= k]

ParetoDistribution/: Domain[ParetoDistribution[k_, alpha_]] :=
	 Interval[{k, Infinity}] /; ParameterQ[ParetoDistribution[k, alpha]]

ParetoDistribution/: PDF[ParetoDistribution[k_, alpha_], x_?Negative] := 0 /;
	ParameterQ[ParetoDistribution[k, alpha]]
ParetoDistribution/: PDF[ParetoDistribution[k_, alpha_], x_] :=
        alpha k^alpha x^(-alpha-1) /; ParameterQ[ParetoDistribution[k, alpha]]

ParetoDistribution/: CDF[ParetoDistribution[k_, alpha_], x_?Negative] := 0 /;
	ParameterQ[ParetoDistribution[k, alpha]]
ParetoDistribution/: CDF[ParetoDistribution[k_, alpha_], x_] :=
        1 - (k/x)^alpha /; ParameterQ[ParetoDistribution[k, alpha]]

ParetoDistribution/: Mean[ParetoDistribution[k_, alpha_]] :=
        alpha k / (alpha - 1) /; ParameterQ[ParetoDistribution[k, alpha]] &&
							 !TrueQ[alpha <= 1]

ParetoDistribution/: Variance[ParetoDistribution[k_, alpha_]] :=
        alpha k^2 / (alpha - 1)^2 / (alpha - 2) /; ParameterQ[
		ParetoDistribution[k, alpha]] && !TrueQ[alpha <= 2]

ParetoDistribution/: StandardDeviation[ParetoDistribution[k_, alpha_]] :=
	Sqrt[alpha k^2 / (alpha - 1)^2 / (alpha - 2)] /; ParameterQ[
		ParetoDistribution[k, alpha]] && !TrueQ[alpha <= 2]

ParetoDistribution/: Skewness[ParetoDistribution[k_, alpha_]] :=
	Module[{mean = alpha k / (alpha - 1),
	        var = alpha k^2 / (alpha - 1)^2 / (alpha - 2), mu3},
	    mu3 = alpha k^3/(alpha - 3) - 3 mean var - mean^3;
	    mu3 / var^(3/2)  
	] /; ParameterQ[ParetoDistribution[k, alpha]] && !TrueQ[alpha <= 3]
 
ParetoDistribution/: Kurtosis[ParetoDistribution[k_, alpha_]] :=
	Module[{mean = alpha k / (alpha - 1),
		var = alpha k^2 / (alpha - 1)^2 / (alpha - 2), mu3, mu4},
	    mu3 = alpha k^3/(alpha - 3) - 3 mean var - mean^3;
	    mu4 = alpha k^4/(alpha - 4) - 4 mean mu3 - 6 mean^2 var - mean^4;
	    mu4 / var^2 
	] /; ParameterQ[ParetoDistribution[k, alpha]] && !TrueQ[alpha <= 4]

ParetoDistribution/: Kurtosis[ParetoDistribution[k_, alpha_]] :=
	Module[{mean = alpha k / (alpha - 1),
		var = alpha k^2 / (alpha - 1)^2 / (alpha - 2), mu3, mu4},
	    mu3 = alpha k^3/(alpha - 3) - 3 mean var - mean^3;
	    mu4 = alpha k^4/(alpha - 4) - 4 mean mu3 - 6 mean^2 var - mean^4;
	    mu4 / var^2 - 3 
	] /; ParameterQ[ParetoDistribution[k, alpha]] && !TrueQ[alpha <= 4]

ParetoDistribution/: CharacteristicFunction[
        ParetoDistribution[k_, alpha_], t_] :=
  (
  (2^alpha*alpha*k^alpha*(t^2)^(alpha/2)*Cos[(alpha*Pi)/2]*Gamma[-alpha])/
    4^(alpha/2) + HypergeometricPFQ[{-alpha/2}, {1/2, 1 - alpha/2}, 
    -(k^2*t^2)/4] - (I*alpha*k*(t^2)^(1/2)*
      HypergeometricPFQ[{1/2 - alpha/2}, {3/2, 3/2 - alpha/2}, 
       -(k^2*t^2)/4]*Sign[t])/(1 - alpha) + 
   (I*2^alpha*k^alpha*(t^2)^(alpha/2)*Gamma[1 - alpha]*Sign[t]*
      Sin[(alpha*Pi)/2])/4^(alpha/2)
  ) /; ParameterQ[ParetoDistribution[k, alpha]]

   (* Integrate[PDF[ParetoDistribution[k, alpha], x] Exp[I x t],
        {x, k, Infinity}, Assumptions -> {Im[t]==0, k > 0, alpha > 0}] *)
 
ParetoDistribution/: Quantile[ParetoDistribution[k_, alpha_], q_] :=
   (
	k / (1-q)^(1/alpha)
   ) /; QuantileQ[q] && ParameterQ[ParetoDistribution[k, alpha]]

ParetoDistribution/: ExpectedValue[f_Function, ParetoDistribution[k_, alpha_],
	opts___?OptionQ] :=
  Module[{x, integral,
	  assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	 integral = Integrate[ f[x] PDF[ParetoDistribution[k, alpha], x],
				{x, k, Infinity},
				Assumptions -> Join[{k>0, alpha>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[ParetoDistribution[k, alpha]]
	
ParetoDistribution/: ExpectedValue[f_, ParetoDistribution[k_, alpha_], x_Symbol,
	opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	 integral = Integrate[ f PDF[ParetoDistribution[k, alpha], x],
		{x, k, Infinity}, Assumptions -> Join[{k>0, alpha>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[ParetoDistribution[k, alpha]]
				
(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE: do not condition Random on ParameterQ so that Random will always
        evaluate and be efficient. *)

pareto = Compile[{{k, _Real}, {alpha, _Real}, {q, _Real}},
		k / (1-q)^(1/alpha)	]	
 
ParetoDistribution/: Random[ParetoDistribution[k_, alpha_]] :=
	pareto[k, alpha, Random[]]
 
(* NOTE: For ParetoDistribution, RandomArray provides no speedup over Random. *)
ParetoDistribution/: RandomArray[ParetoDistribution[k_, alpha_], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[pareto[k, alpha, Random[]], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* ========================== Uniform Distribution ======================= *)
(* Uniform Distribution, K & J, Vol. 2, Chap. 25 *)

UniformDistribution/: ParameterQ[UniformDistribution[min_, max_]] := And[
        If[FreeQ[N[max], Complex], True,
           Message[UniformDistribution::realparm, max]; False],
        If[FreeQ[N[min], Complex], True,
           Message[UniformDistribution::realparm, min]; False],
        If[N[max-min] > 0, True,
           Message[UniformDistribution::posparm, max-min]; False, True]
]

UniformDistribution::posparm =
"The parameter difference (max-min) = `` is expected to be positive."

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

UniformDistribution/: DomainQ[UniformDistribution[min_, max_], list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
	Scan[If[!TrueQ[min <= # <= max], Return[False]]&, list] =!= False
UniformDistribution/: DomainQ[UniformDistribution[min_, max_], x_] :=
		FreeQ[N[x], Complex] && TrueQ[min <= x <= max]

UniformDistribution/: Domain[UniformDistribution[min_:0, max_:1]] :=
        Interval[{min, max}] /; ParameterQ[UniformDistribution[min, max]]

UniformDistribution/: PDF[UniformDistribution[min_:0, max_:1], x_] :=
    (Sign[x-min] - Sign[x-max])/(2 (max - min)) /; ParameterQ[
					UniformDistribution[min, max]]

UniformDistribution/: CDF[UniformDistribution[min_:0, max_:1], x_] :=
    (1+Sign[x-min] (x-min)/(max-min)-(max-x) Sign[max-x]/
	(max-min))/2 		/; ParameterQ[UniformDistribution[min, max]]

UniformDistribution/: Mean[UniformDistribution[min_:0, max_:1]] := 
	(max+min)/2		/; ParameterQ[UniformDistribution[min, max]]

UniformDistribution/: Variance[UniformDistribution[min_:0, max_:1]] :=
	(max-min)^2 / 12	/; ParameterQ[UniformDistribution[min, max]]

UniformDistribution/: StandardDeviation[UniformDistribution[min_:0, max_:1]] :=
    (max-min)/Sqrt[12]		/; ParameterQ[UniformDistribution[min, max]]

UniformDistribution/: Skewness[UniformDistribution[___]] = 0
 
UniformDistribution/: Kurtosis[UniformDistribution[___]] = 9/5 

UniformDistribution/: KurtosisExcess[UniformDistribution[___]] = -6/5 

UniformDistribution/: CharacteristicFunction[
                UniformDistribution[min_:0, max_:1], t_] :=
        (Exp[I max t] - Exp[I min t])/((max - min) I t) /; ParameterQ[
						UniformDistribution[min, max]]

UniformDistribution/: Quantile[UniformDistribution[min_:0, max_:1], q_] :=
  q max + (1-q) min /; QuantileQ[q] && ParameterQ[UniformDistribution[min, max]]

UniformDistribution/: ExpectedValue[f_Function,
			UniformDistribution[min_:0, max_:1], opts___?OptionQ] :=
 Module[{x, integral,
	assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f[x], {x, min, max},
			Assumptions -> assmp]/(max-min);
	 FreeQ[integral, Integrate])
 ] /; ParameterQ[UniformDistribution[min, max]]

UniformDistribution/: ExpectedValue[f_, UniformDistribution[min_:0, max_:1],
	x_Symbol, opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (integral = Integrate[ f, {x, min, max},
			Assumptions -> assmp]/(max-min);
         FreeQ[integral, Integrate])
  ] /; ParameterQ[UniformDistribution[min, max]]

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

UniformDistribution/: Random[UniformDistribution[min_:0, max_:1]] :=
		Random[Real, {min, max}]

(* NOTE: For UniformDistribution, RandomArray provides no speedup over Random. *)
UniformDistribution/: RandomArray[UniformDistribution[min_:0, max_:1], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[Random[Real, {min, max}], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]
 

(* =========================== Weibull Distribution ====================== *)
(* Weibull Distribution, K & J, Vol. 1, Chap. 20 *)

WeibullDistribution/: ParameterQ[WeibullDistribution[alpha_, beta_]] := And[
        If[FreeQ[N[alpha], Complex], True,
           Message[WeibullDistribution::realparm, alpha]; False],
        If[FreeQ[N[beta], Complex], True,
           Message[WeibullDistribution::realparm, beta]; False],
        If[N[alpha] > 0, True,
           Message[WeibullDistribution::posparm, alpha]; False, True],
        If[N[beta] > 0, True,
           Message[WeibullDistribution::posparm, beta]; False, True]
]

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

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

WeibullDistribution/: DomainQ[WeibullDistribution[alpha_, beta_],
	 list_?VectorQ] :=
	FreeQ[N[list], Complex] &&
        Scan[If[!TrueQ[# >= 0], Return[False]]&, list] =!= False
WeibullDistribution/: DomainQ[WeibullDistribution[alpha_, beta_], x_] :=
	FreeQ[N[x], Complex] && TrueQ[x >= 0]

WeibullDistribution/: Domain[WeibullDistribution[___]] := 
		Interval[{0, Infinity}]

WeibullDistribution/: PDF[WeibullDistribution[alpha_, beta_], x_?Negative] :=
				0 /; ParameterQ[WeibullDistribution[alpha, beta]]
WeibullDistribution/: PDF[WeibullDistribution[alpha_, beta_], x_] :=
        alpha x^(alpha-1) Exp[-(x/beta)^alpha] / beta^alpha /; ParameterQ[
				WeibullDistribution[alpha, beta]]

WeibullDistribution/: CDF[WeibullDistribution[alpha_, beta_], x_?Negative] :=
				0 /; ParameterQ[WeibullDistribution[alpha, beta]]
WeibullDistribution/: CDF[WeibullDistribution[alpha_, beta_], x_] :=
        1 - Exp[-(x/beta)^alpha] /; ParameterQ[WeibullDistribution[alpha, beta]]

WeibullDistribution/: Mean[WeibullDistribution[alpha_, beta_]] :=
	beta Gamma[1+1/alpha] /; ParameterQ[WeibullDistribution[alpha, beta]]

WeibullDistribution/: Variance[WeibullDistribution[alpha_, beta_]] :=
	beta^2 (Gamma[1+2/alpha] - Gamma[1+1/alpha]^2) /; ParameterQ[
					WeibullDistribution[alpha, beta]]

WeibullDistribution/: StandardDeviation[
                WeibullDistribution[alpha_, beta_]] :=
	beta Sqrt[Gamma[1+2/alpha] - Gamma[1+1/alpha]^2] /; ParameterQ[
					WeibullDistribution[alpha, beta]]

WeibullDistribution/: Skewness[WeibullDistribution[alpha_, beta_]] :=
        (Gamma[1+3/alpha] - 3 Gamma[1+1/alpha] Gamma[1+2/alpha] +
            2 Gamma[1+1/alpha]^3 ) /
            (Gamma[1+2/alpha] - Gamma[1+1/alpha]^2)^(3/2) /; ParameterQ[
					WeibullDistribution[alpha, beta]]

WeibullDistribution/: Kurtosis[WeibullDistribution[alpha_, beta_]] :=
        (Gamma[1+4/alpha] - 4 Gamma[1+1/alpha] Gamma[1+3/alpha] +
            6 Gamma[1+1/alpha]^2 Gamma[1+2/alpha] -
            3 Gamma[1+1/alpha]^4) /
            (Gamma[1+2/alpha] - Gamma[1+1/alpha]^2)^2 /; ParameterQ[
					WeibullDistribution[alpha, beta]]

WeibullDistribution/: KurtosisExcess[WeibullDistribution[alpha_, beta_]] :=
        (Gamma[1+4/alpha] - 4 Gamma[1+1/alpha] Gamma[1+3/alpha] +
            6 Gamma[1+1/alpha]^2 Gamma[1+2/alpha] -
            3 Gamma[1+1/alpha]^4) /
            (Gamma[1+2/alpha] - Gamma[1+1/alpha]^2)^2 - 3 /; ParameterQ[
					WeibullDistribution[alpha, beta]]

WeibullDistribution/: CharacteristicFunction[
        WeibullDistribution[alpha_, beta_], t_] :=
   Integrate[PDF[WeibullDistribution[alpha, beta], x] Exp[I x t], 
	{x, 0, Infinity}, Assumptions -> {Im[t]==0}]	/; False
 
WeibullDistribution/: Quantile[WeibullDistribution[alpha_, beta_], q_] :=
	beta (-Log[1 - q])^(1/alpha) /; QuantileQ[q] && ParameterQ[
					WeibullDistribution[alpha, beta]]

WeibullDistribution/: ExpectedValue[f_Function,
			 WeibullDistribution[alpha_, beta_], opts___?OptionQ] :=
  Module[{x, integral,
	  assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	integral = Integrate[ f[x] PDF[WeibullDistribution[alpha, beta], x],
			{x, 0, Infinity},
			Assumptions -> Join[{alpha>0, beta>0}, assmp]];
	 FreeQ[integral, Integrate])
  ] /; ParameterQ[WeibullDistribution[alpha, beta]]

WeibullDistribution/: ExpectedValue[f_, WeibullDistribution[alpha_, beta_],
		x_Symbol, opts___?OptionQ] :=
  Module[{integral, assmp = Assumptions /. {opts} /. Options[ExpectedValue]},
   (
	integral
   ) /; (
	integral = Integrate[ f PDF[WeibullDistribution[alpha, beta], x],
                        {x, 0, Infinity},
                        Assumptions -> Join[{alpha>0, beta>0}, assmp]];
         FreeQ[integral, Integrate])
  ] /; ParameterQ[WeibullDistribution[alpha, beta]]

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

weibull = Compile[{{alpha, _Real}, {beta, _Real}, {q, _Real}},
		beta Log[1/q]^(1/alpha)		]

WeibullDistribution/: Random[WeibullDistribution[alpha_, beta_]] :=
	weibull[alpha, beta, Random[]]

(* NOTE: For WeibullDistribution, RandomArray provides no speedup over Random. *)
WeibullDistribution/: RandomArray[WeibullDistribution[alpha_, beta_], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Table[weibull[alpha, beta, Random[]], {n}];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],
       array  ]
  ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]
 
		

End[]

SetAttributes[ BetaDistribution, ReadProtected];
SetAttributes[ CauchyDistribution, ReadProtected];
SetAttributes[ ChiDistribution, ReadProtected];
SetAttributes[ NoncentralChiSquareDistribution, ReadProtected];
SetAttributes[ ExponentialDistribution, ReadProtected];
SetAttributes[ ExtremeValueDistribution, ReadProtected];
SetAttributes[ NoncentralFRatioDistribution, ReadProtected];
SetAttributes[ GammaDistribution, ReadProtected];
SetAttributes[ HalfNormalDistribution, ReadProtected];
SetAttributes[ LaplaceDistribution, ReadProtected];
SetAttributes[ LogNormalDistribution, ReadProtected];
SetAttributes[ LogisticDistribution, ReadProtected];
SetAttributes[ RayleighDistribution, ReadProtected];
SetAttributes[ NoncentralStudentTDistribution, ReadProtected];
SetAttributes[ ParetoDistribution, ReadProtected];
SetAttributes[ UniformDistribution, ReadProtected];
SetAttributes[ WeibullDistribution, ReadProtected];


Protect[BetaDistribution,CauchyDistribution,ChiDistribution,
	NoncentralChiSquareDistribution,ExponentialDistribution, 
	ExtremeValueDistribution,NoncentralFRatioDistribution,
	GammaDistribution,HalfNormalDistribution, LaplaceDistribution,
	LogNormalDistribution,LogisticDistribution,
	RayleighDistribution,NoncentralStudentTDistribution, 
	ParetoDistribution, UniformDistribution, WeibullDistribution];

EndPackage[]

