(*:Mathematica Version: 3.0 *)

(*:Package Version: 1.1 *)

(*:Name: Statistics`NormalDistribution` *)

(*:Context: Statistics`NormalDistribution` *)

(*:Title: Statistical Distributions Derived from the Normal Distribution *)

(*:Author:
  David Withoff (Wolfram Research), March 1991
*)

(*:History:
  David Withoff (Wolfram Research), March 1991.
  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, and RandomArray,
	 speeded up random number generation, redefined Domain in terms
	 of Interval, ECM (Wolfram Research), February 1995.
*)

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

(*:Reference: Usage messages only. *)

(*:Summary:
This package provides properties and functionals of the four standard
probability distributions derived from the normal (Gaussian) distribution.
The four distributions are the normal distribution, Student's t-distribution,
the chi-square distribution, and the F-ratio distribution.
*)

(*:Keywords: Continuous distribution, normal 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.
*)

(*:Note: Most functions provide numerical results. *)
 
(*: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`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."

]


        (*  Distributions  *)

ChiSquareDistribution::usage =
"ChiSquareDistribution[n] represents the chi-square distribution
with n degrees of freedom."

FRatioDistribution::usage =
"FRatioDistribution[n1, n2] represents the F-ratio distribution with n1
numerator degrees of freedom and n2 denominator degrees of freedom."

NormalDistribution::usage =
"NormalDistribution[mu, sigma] represents the normal (Gaussian)
distribution with mean mu and standard deviation sigma."

StudentTDistribution::usage =
"StudentTDistribution[n] represents Student's T distribution
with n degrees of freedom."

(* Introduce usage messages for new functions *)
 
PercentagePoint::usage =
"PercentagePoint[distribution, q] gives the percentage point corresponding to
percentage q of the specified statistical distribution.  This is equivalent to
Quantile[distribution, (q+1)/2] for symmetrical distributions."

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

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

ExpectedValue::usage = ExpectedValue::usage <> " " <>
"ExpectedValue[f, distribution] gives the expected value of the pure function
f with respect to the specified statistical distribution.
ExpectedValue[f, distribution, x] gives the expected value of the
function f of x with respect to the specified univariate statistical
distribution for x.";

Mean::usage = Mean::usage <> " " <>
"Mean[distribution] gives the mean of the specified statistical distribution.";

Variance::usage = Variance::usage <> " " <>
"Variance[distribution] gives the variance of the specified statistical
distribution.";

StandardDeviation::usage = StandardDeviation::usage <> " " <>
"StandardDeviation[distribution] gives the standard deviation of the
specified statistical distribution.";

Skewness::usage = Skewness::usage <> " " <>
"Skewness[distribution] gives the coefficient of skewness of the
specified statistical distribution.";

Kurtosis::usage = Kurtosis::usage <> " " <>
"Kurtosis[distribution] gives the coefficient of kurtosis of the
specified statistical distribution.";

KurtosisExcess::usage = KurtosisExcess::usage <> " " <>
"KurtosisExcess[distribution] gives the kurtosis excess for the
specified statistical distribution.";

Quantile::usage = Quantile::usage <> " " <>
"Quantile[distribution, q] gives the q-th quantile of the specified
statistical distribution."

]
 
Unprotect[ChiSquareDistribution, FRatioDistribution,
        NormalDistribution, StudentTDistribution];

Begin["`Private`"]

(* users may set Assumptions option of ExpectedValue, which is then passed
	on to Integrate, in addition to the built-in assumptions
	regarding parameters *)
Options[ExpectedValue] = {Assumptions -> {}}

(* ========================= Chi-square Distribution ======================= *)
(* J & K, CUD Vol. 1, Chap. 17 *)

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

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

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

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

ChiSquareDistribution/: Domain[ChiSquareDistribution[n_]] :=
	Interval[{0, Infinity}] /; ParameterQ[ChiSquareDistribution[n]]

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

ChiSquareDistribution/: Mean[ChiSquareDistribution[n_]] := 
	n /; ParameterQ[ChiSquareDistribution[n]]

ChiSquareDistribution/: Variance[ChiSquareDistribution[n_]] := 
	2 n /; ParameterQ[ChiSquareDistribution[n]]

ChiSquareDistribution/: StandardDeviation[ChiSquareDistribution[n_]] :=
	Sqrt[2 n] /; ParameterQ[ChiSquareDistribution[n]]

ChiSquareDistribution/: Skewness[ChiSquareDistribution[n_]] := 
	2 Sqrt[2/n] /; ParameterQ[ChiSquareDistribution[n]]

ChiSquareDistribution/: Kurtosis[ChiSquareDistribution[n_]] := 
	3 + 12/n /; ParameterQ[ChiSquareDistribution[n]]

ChiSquareDistribution/: KurtosisExcess[ChiSquareDistribution[n_]] :=
	12/n /; ParameterQ[ChiSquareDistribution[n]]

ChiSquareDistribution/: CharacteristicFunction[
                ChiSquareDistribution[n_], t_] :=
        (1 - 2 I t)^(-n/2) /; ParameterQ[ChiSquareDistribution[n]]

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

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

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

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

ChiSquareDistribution/: Random[ChiSquareDistribution[n_]] :=
  Which[n < 2,
          2 gsGamma[n/2, n/2/E],
	n == 2,
	  2 exponential[], 
        n < 5,
          a = n/2-1;
          2 a gkm1Gamma[(n/2 - 1/(3 n))/a, 2/a],
        True,
          a = n/2-1;
          2 a gkm2Gamma[1/Sqrt[n/2], (n/2 - 1/(3 n))/a, 2/a]
  ]

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

ChiSquareDistribution/: RandomArray[ChiSquareDistribution[n_], dim_] :=
  Module[{m, array, r, k, a, b, c},
    m = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Which[n < 2,
                      r = n/2/E;
                      2 Table[gsGamma[n/2, r], {m}],
		  n == 2,
		      2 Table[exponential[], {m}],
                  n < 5,
                      a = n/2-1;  b = (n/2 - 1/(6 n/2))/a;  c = 2/a;
                      2 a Table[gkm1Gamma[b, c], {m}],
                  True,
                      k = 1/Sqrt[n/2];
                      a = n/2-1;  b = (n/2 - 1/(6 n/2))/a;  c = 2/a;
                      2 a Table[gkm2Gamma[k, b, c], {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 = Compile[{}, -Log[Random[]]]

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* 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. *)
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. *)
gkm1Gamma = Compile[{{b, _Real}, {c, _Real}},
   Module[{u1 = Random[], u2 = Random[], w},
     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. *)
gkm2Gamma = Compile[{{k, _Real}, {b, _Real}, {c, _Real}},
   Module[{u1 = Random[], u2, w},
     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
   ]    ]


(* Suppose x is distributed according to ChiSquareDistribution[n].
 Then as n -> Infinity, the distribution of (x-Mean[x])/StandardDeviation[x]
	approaches the NormalDistribution[0, 1] distribution.
 How best to represent this using ChiSquareDistribution[Infinity] ?
*)

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

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

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

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

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

FRatioDistribution/: Domain[FRatioDistribution[n1_, n2_]] :=
	Interval[{0, Infinity}] /; ParameterQ[FRatioDistribution[n1, n2]]

FRatioDistribution/: PDF[FRatioDistribution[n1_, n2_], x_?Negative] := 0 /;
				ParameterQ[FRatioDistribution[n1, n2]]
FRatioDistribution/: PDF[FRatioDistribution[n1_, n2_], x_] :=
	n1^(n1/2) n2^(n2/2) x^(n1/2 - 1)/
            		((n2 + n1 x)^((n1 + n2)/2) Beta[n1/2, n2/2]) /;
				ParameterQ[FRatioDistribution[n1, n2]]

FRatioDistribution/: CDF[FRatioDistribution[n1_, n2_], x_?Negative] := 0 /;
				ParameterQ[FRatioDistribution[n1, n2]]
FRatioDistribution/: CDF[FRatioDistribution[n1_, n2_], x_] :=
	BetaRegularized[n2/(n2 + n1 x), 1, n2/2, n1/2] /; ParameterQ[
						FRatioDistribution[n1, n2]]

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

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

FRatioDistribution/: StandardDeviation[FRatioDistribution[n1_, n2_]] :=
   (n2/(n2-2)) Sqrt[2] Sqrt[n1+n2-2]/(Sqrt[n1] Sqrt[n2-4]) /;
		ParameterQ[FRatioDistribution[n1, n2]] && !TrueQ[n2 <= 4]

FRatioDistribution/: Skewness[FRatioDistribution[n1_, n2_]] :=
	(2 n1 + n2 - 2) Sqrt[8 (n2-4)] / (Sqrt[n1] (n2-6) Sqrt[n1+n2-2]) /;
		ParameterQ[FRatioDistribution[n1, n2]] && !TrueQ[n2 <= 6]

FRatioDistribution/: Kurtosis[FRatioDistribution[n1_, n2_]] :=
	12 ((n2-2)^2 (n2-4) + n1 (n1+n2-2) (5 n2 -22))/
		(n1 (n2-6) (n2-8) (n1+n2-2)) + 3 /;
		ParameterQ[FRatioDistribution[n1, n2]] && !TrueQ[n2 <= 8]

FRatioDistribution/: KurtosisExcess[FRatioDistribution[n1_, n2_]] :=
	12 ((n2-2)^2 (n2-4) + n1 (n1+n2-2) (5 n2 -22))/
	   (n1 (n2-6) (n2-8) (n1+n2-2)) /; ParameterQ[
			FRatioDistribution[n1, n2]] && !TrueQ[n2 <= 8]

FRatioDistribution/:
        CharacteristicFunction[FRatioDistribution[n1_, n2_], t_] :=
    Hypergeometric1F1[n1/2, 1 - n2/2, -I t n2/n1] /; ParameterQ[
						FRatioDistribution[n1, n2]]

FRatioDistribution/: Quantile[FRatioDistribution[n1_, n2_], q_] :=
    n2/n1 (1/InverseBetaRegularized[1, -q,
		 n2/2, n1/2] - 1) /; QuantileQ[q] && ParameterQ[
					FRatioDistribution[n1, n2]]

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

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

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

FRatioDistribution/: Random[FRatioDistribution[n1_, n2_]] :=
    n2/n1 Random[ChiSquareDistribution[n1]]/Random[ChiSquareDistribution[n2]]

FRatioDistribution/: RandomArray[FRatioDistribution[n1_, n2_], dim_] :=
   (
    n2/n1 RandomArray[ChiSquareDistribution[n1], dim] /
	  RandomArray[ChiSquareDistribution[n2], dim]
   ) /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]


(* Suppose f is distributed according to FRatioDistribution[n1, n2].
	As n1 -> Infinity, the distribution of n2/f approaches the
		ChiSquareDistribution[n2] distribution.
	As n2 -> Infinity, the distribution of n1*f approaches the
		ChiSquareDistribution[n1] distribution. 
 How best to represent this using FRatioDistribution[Infinity, n2] or
	FRatioDistribution[n1, Infinity] ?
*)


(* ========================= Normal Distribution ========================= *)
(* K & J, CUD Vol. 1, Chap. 13 *)

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

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

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

NormalDistribution/: DomainQ[NormalDistribution[mu_, sigma_], x_] :=
	 FreeQ[N[x], Complex]

NormalDistribution/: Domain[NormalDistribution[mu_:0, sigma_:1]] :=
  Interval[{-Infinity, Infinity}] /; ParameterQ[NormalDistribution[mu, sigma]]

NormalDistribution/: PDF[NormalDistribution[mu_:0, sigma_:1], x_] :=
        Exp[-((x-mu)/sigma)^2 / 2]/(sigma Sqrt[2Pi]) /; ParameterQ[
				NormalDistribution[mu, sigma]]

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

NormalDistribution/: Mean[NormalDistribution[mu_:0, sigma_:1]] := 
	mu  /; ParameterQ[NormalDistribution[mu, sigma]]

NormalDistribution/: StandardDeviation[
                NormalDistribution[mu_:0, sigma_:1]] := 
	sigma /; ParameterQ[NormalDistribution[mu, sigma]]

NormalDistribution/: Variance[NormalDistribution[mu_:0, sigma_:1]] := 
	sigma^2 /; ParameterQ[NormalDistribution[mu, sigma]]

NormalDistribution/: Skewness[NormalDistribution[mu_:0, sigma_:1]] = 0

NormalDistribution/: Kurtosis[NormalDistribution[mu_:0, sigma_:1]] = 3
 
NormalDistribution/: KurtosisExcess[NormalDistribution[mu_:0, sigma_:1]] = 0
 
NormalDistribution/: CharacteristicFunction[
                NormalDistribution[mu_:0, sigma_:1], t_] :=
      	Exp[I mu t - sigma^2 t^2/2] /; ParameterQ[
					NormalDistribution[mu, sigma]]
 
NormalDistribution/: Quantile[NormalDistribution[mu_:0, sigma_:1], q_] :=
	mu + Sqrt[2] sigma InverseErf[2 q - 1] /; QuantileQ[q] &&
				ParameterQ[NormalDistribution[mu, sigma]]

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

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

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* NOTE: do not condition Random on ParameterQ so that Random will always
        evaluate and be efficient. *)
 
normal = Compile[{{mu, _Real}, {sigma, _Real}, {q1, _Real}, {q2, _Real}},
		mu + sigma Sqrt[-2 Log[q1]] Cos[2Pi q2]	]

NormalDistribution/: Random[NormalDistribution[mu_:0, sigma_:1]] :=
	normal[mu, sigma, Random[], Random[]]

(*
Now we do not need to compile this,  the
Table command will auto-compile.

normalpair = Compile[{{mu, _Real}, {sigma, _Real}, {q1, _Real}, {q2, _Real}},
		mu + sigma Sqrt[-2 Log[q1]] {Cos[2Pi q2], Sin[2Pi q2]}   ]
*)


NormalDistribution/: RandomArray[NormalDistribution[mu_:0, sigma_:1], dim_] :=
  Module[{n, array},
    n = If[VectorQ[dim], Apply[Times, dim], dim];
    array = Flatten[
    	Table[
    	    Block[{q1 = Random[], q2 = Random[]},
    	    	mu + sigma Sqrt[-2. Log[q1]] {Cos[2. Pi q2], 
		Sin[ 2. Pi q2]}], {Quotient[ n, 2]}]];
    If[OddQ[n],
       AppendTo[array, normal[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)&]


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

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

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

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

StudentTDistribution/: DomainQ[StudentTDistribution[n_], x_] :=
	 FreeQ[N[x], Complex]

StudentTDistribution/: Domain[StudentTDistribution[n_]] :=
        Interval[{-Infinity, Infinity}] /; ParameterQ[StudentTDistribution[n]]

StudentTDistribution/: PDF[StudentTDistribution[n_], x_] :=
      1/(Sqrt[n] Beta[n/2, 1/2]) Sqrt[n/(n+x^2)]^(n+1) /;
				ParameterQ[StudentTDistribution[n]]

StudentTDistribution/: CDF[StudentTDistribution[n_], x_] := 
      (1 + Sign[x] BetaRegularized[n/(n+x^2), 1, n/2, 1/2])/2 /;
				ParameterQ[StudentTDistribution[n]]

StudentTDistribution/: Mean[StudentTDistribution[n_]] :=
	Indeterminate /; ParameterQ[StudentTDistribution[n]] && TrueQ[n <= 1]

StudentTDistribution/: Mean[StudentTDistribution[n_]] :=
	0 /; ParameterQ[StudentTDistribution[n]] && !TrueQ[n <= 1]

StudentTDistribution/: Variance[StudentTDistribution[n_]] :=
	Indeterminate /; ParameterQ[StudentTDistribution[n]] && TrueQ[n <= 2]

StudentTDistribution/: Variance[StudentTDistribution[n_]] :=
	n/(n-2) /; ParameterQ[StudentTDistribution[n]] && !TrueQ[n <= 2]

StudentTDistribution/: StandardDeviation[StudentTDistribution[n_]] :=
	Indeterminate /; ParameterQ[StudentTDistribution[n]] && TrueQ[n <= 2]

StudentTDistribution/: StandardDeviation[StudentTDistribution[n_]] :=
	Sqrt[n/(n-2)] /; ParameterQ[StudentTDistribution[n]] && !TrueQ[n <= 2]

StudentTDistribution/: Skewness[StudentTDistribution[n_]] :=
	Indeterminate /; ParameterQ[StudentTDistribution[n]] && TrueQ[n <= 3]

StudentTDistribution/: Skewness[StudentTDistribution[n_]] :=
	0 /; ParameterQ[StudentTDistribution[n]] && !TrueQ[n <= 3]

StudentTDistribution/: Kurtosis[StudentTDistribution[n_]] :=
	Indeterminate /; ParameterQ[StudentTDistribution[n]] && TrueQ[n <= 4]

StudentTDistribution/: Kurtosis[StudentTDistribution[n_]] :=
      	6/(n-4) + 3 /; ParameterQ[StudentTDistribution[n]] && !TrueQ[n <= 4]

StudentTDistribution/: KurtosisExcess[StudentTDistribution[n_]] := 
      	Indeterminate /; ParameterQ[StudentTDistribution[n]] && TrueQ[n <= 4]

StudentTDistribution/: KurtosisExcess[StudentTDistribution[n_]] := 
      	6/(n-4) /; ParameterQ[StudentTDistribution[n]] && !TrueQ[n <= 4]

StudentTDistribution/:
        CharacteristicFunction[StudentTDistribution[1], t_] := Exp[-t Sign[t]]

StudentTDistribution/:
        CharacteristicFunction[StudentTDistribution[n_], t_] :=
      2 BesselK[n/2, Abs[t] Sqrt[n]] (Abs[t] Sqrt[n]/2)^(n/2) /
		Gamma[n/2] /; ParameterQ[StudentTDistribution[n]]

StudentTDistribution/: Quantile[StudentTDistribution[n_], q_] :=
      Sign[2 q - 1] Sqrt[1/InverseBetaRegularized[
        1, (1 - 2 q) Sign[2 q - 1], n/2, 1/2] -1] Sqrt[n] /; QuantileQ[q] &&
				ParameterQ[StudentTDistribution[n]]

(* Percentage point table in Abramowitz and Stegun *)
StudentTDistribution/: PercentagePoint[StudentTDistribution[n_], q_] :=
   Sqrt[n] Sqrt[1/InverseBetaRegularized[1, -q, n/2, 1/2] - 1] /; ParameterQ[
				StudentTDistribution[n]]

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

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

StudentTDistribution/: Random[StudentTDistribution[n_]] :=
	Random[NormalDistribution[0, 1]] /
		 Sqrt[Random[ChiSquareDistribution[n]]/n]

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


(* Suppose z is distributed according to StudentTDistribution[n].  As
	n -> Infinity, the distribution of z approaches the
	NormalDistribution[0, 1] distribution.
*)

StudentTDistribution[Infinity] = NormalDistribution[0, 1]

(* ===================================================================== *)
End[]
SetAttributes[ ChiSquareDistribution, ReadProtected];
SetAttributes[ FRatioDistribution, ReadProtected];    
SetAttributes[ NormalDistribution, ReadProtected];    
SetAttributes[ StudentTDistribution, ReadProtected];    

Protect[ChiSquareDistribution, FRatioDistribution,
        NormalDistribution, StudentTDistribution];

EndPackage[]
