(*:Mathematica Version: 2.0 *)

(*:Package Version: 1.2 *)

(*:Name: Statistics`ConfidenceIntervals` *)

(*:Context: Statistics`ConfidenceIntervals` *)

(*:Title: Confidence Intervals Related to the Normal Distribution *)

(*:Author: Dave Withoff *)

(*:History:
  Dave Withoff (Wolfram Research), February 1990.
  Revised in March 1991, July 1995.
*)

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

(*:Reference: Usage messages only. *)

(*:Summary:
This package provides confidence intervals for elementary reference
distributions derived from the normal distribution.  Distributions
represented are NormalDistribution, StudentTDistribution,
ChiSquareDistribution, and FRatioDistribution.
*)

(*:Keywords: Confidence interval, confidence level *)

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

(*:Warning: None. *)

(*:Sources: Basic statistics texts *)

BeginPackage["Statistics`ConfidenceIntervals`",
             "Statistics`DescriptiveStatistics`",
             "Statistics`NormalDistribution`",
	     "Statistics`Common`PopulationsCommon`"]
 
MeanCI::usage =
"MeanCI[list, options] returns a list {lower, upper}
representing a confidence interval for the population mean, using
the entries in list as a sample drawn from the population."

Options[MeanCI] = {ConfidenceLevel -> .95,
    KnownStandardDeviation -> None, KnownVariance -> None}

VarianceCI::usage =
"VarianceCI[list] returns a list {lower, upper} representing
a confidence interval for the population variance, using the entries in
list as a sample drawn from the population."

Options[VarianceCI] = {ConfidenceLevel -> .95}

MeanDifferenceCI::usage =
"MeanDifferenceCI[list1, list2, options] returns a list {lower, upper}
representing a confidence interval for the difference Mean[list1]
- Mean[list2], using the entries in list1 as a sample from the
first population and the entries in list2 as a sample from the
second population."

Options[MeanDifferenceCI] = {ConfidenceLevel -> .95,
    KnownStandardDeviation -> None, KnownVariance -> None,
    EqualVariances -> False}

VarianceRatioCI::usage =
"VarianceRatioCI[list1, list2, options] returns a list representing
a confidence interval for the ratio of estimated population variances,
Variance[list1]/Variance[list2], using the entries in list1 as a
sample from the first population and the entries in list2 as a
sample from the second population."
 
Options[VarianceRatioCI] := {ConfidenceLevel -> .95}

StudentTCI::usage =
"StudentTCI[mean, se, dof, ConfidenceLevel -> c] returns
a list {lower, upper} representing a confidence interval at 
confidence level c for the population mean, based on a sample mean,
its standard error se and dof degrees of freedom. This function is used
by MeanCI when the population variance is estimated from
the sample."

Options[StudentTCI] = {ConfidenceLevel -> .95}

NormalCI::usage =
"NormalCI[mean, sd, ConfidenceLevel -> c] returns a list
{lower, upper} representing a confidence interval at confidence level
c for the population mean, based on the sample mean and its standard
deviation. This function is used by MeanCI when the population 
variance is specified."

Options[NormalCI] = {ConfidenceLevel -> .95}

ChiSquareCI::usage =
"ChiSquareCI[var, dof, ConfidenceLevel -> c] returns a list
{lower, upper} representing a confidence interval at confidence level c
for the population variance, based on a sample with dof degrees of
freedom and unbiased variance estimate var.  This function is used
by VarianceCI."

Options[ChiSquareCI] = {ConfidenceLevel -> .95}

FRatioCI::usage =
"FRatioCI[ratio, numdof, dendof, ConfidenceLevel -> c] returns
a list {lower, upper} representing a confidence interval at confidence
level c for ratio, where ratio is a ratio of sample variances, numdof is the
number of numerator degrees of freedom, and dendof is the number of
denominator degrees of freedom.  This function is used by
VarianceRatioCI."

Options[FRatioCI] = {ConfidenceLevel -> .95}


	(*  Options  *)

ConfidenceLevel::usage =
"ConfidenceLevel is an option used to specify the confidence level
in functions that return statistical confidence intervals."

Begin["`Private`"]

(* Confidence Intervals for one Mean *)

MeanCI[list_List, options___] :=
    Block[{mean = Mean[list], n = Length[list], c, var0, sd0},
        {c, var0, sd0} =
            {ConfidenceLevel, KnownVariance, KnownStandardDeviation} /.
                {options} /. Options[MeanCI];
        If[!(0 <= N[c] <= 1), Message[MeanCI::clev, c]];
        If[sd0 =!= None,
            If[var0 =!= None, Message[MeanCI::varsd]];
            var0 = sd0^2
        ];
        If[var0 === None,
            Return[StudentTCI[mean, StandardErrorOfSampleMean[list],
                                    n-1, ConfidenceLevel->c]],
            Return[NormalCI[mean, Sqrt[var0/n], ConfidenceLevel->c]]
        ]
    ]

MeanCI::varsd = "Warning: KnownStandardDeviation and KnownVariance
have both been specified.  KnownVariance will be ignored."

MeanCI::clev = "Warning: ConfidenceLevel specification `1` is not
in the expected range 0 <= ConfidenceLevel <= 1."

StudentTCI[mean_, se_, dof_, options___] :=
    Block[{c = ConfidenceLevel /. {options} /. Options[StudentTCI], h},
        If[N[se] <= 0, Message[StudentTCI::badse, se]];
        h = se Quantile[StudentTDistribution[dof], (1+c)/2];
        {mean-h, mean+h}]

StudentTCI::badse = "Warning: Standard error `1` is normally
a positive quantity."

NormalCI[mean_, sd_, options___] :=
    Block[{c = ConfidenceLevel /. {options} /. Options[NormalCI], h},
        h = sd Quantile[NormalDistribution[0, 1], (1+c)/2];
        {mean-h, mean+h}]

(* Confidence Intervals for one Variance *) 

VarianceCI[list_List, options___] :=
    Block[{dof = Length[list] - 1, var = Variance[list],
	   c = ConfidenceLevel /. {options} /. Options[VarianceCI]},
        ChiSquareCI[var, dof, ConfidenceLevel->c]
    ]

ChiSquareCI[var_, dof_, options___] :=
    Block[{c = ConfidenceLevel /. {options} /. Options[ChiSquareCI],
           sumsq = dof var, lower, upper},
        upper = sumsq/Quantile[ChiSquareDistribution[dof], (1-c)/2];
        lower = sumsq/Quantile[ChiSquareDistribution[dof], (1+c)/2];
        {lower, upper}
    ]

(* Comparison of two Means *)

MeanDifferenceCI[list1_List, list2_List, options___] :=
    Block[{diff = Mean[list1] -  Mean[list2], var0, sd0},
        {var0, sd0, equalvar, c} = {KnownVariance, KnownStandardDeviation,
	 	EqualVariances, ConfidenceLevel} /.
            {options} /. Options[MeanDifferenceCI];
        If[sd0 =!= None,
            If[var0 =!= None, Message[MeanDifferenceCI::varsd]];
            var0 = sd0^2
        ];
        mdInterval[list1, list2, diff, var0, equalvar, c]
    ]

MeanDifferenceCI::varsd = "Warning: KnownStandardDeviation
and KnownVariance have both been specified.  KnownVariance will be ignored."

mdInterval[list1_, list2_, diff_, None, equalvar_, c_] :=
    Block[{var1 = Variance[list1], var2 = Variance[list2],
	   n1 = Length[list1], n2 = Length[list2], dof, stderr},
        If[TrueQ[equalvar],
            dof = n1 + n2 - 2;
            stderr = Sqrt[(((n1-1) var1 + (n2-1) var2) / dof)
                              (1/n1 + 1/n2)],
        (* else *)
            stderr = Sqrt[var1/n1 + var2/n2];
            dof = stderr^4 /
                  ((var1/n1)^2/(n1-1) + (var2/n2)^2/(n2-1))
        ];
        StudentTCI[diff, stderr, dof, ConfidenceLevel->c]
    ]

mdInterval[list1_, list2_, diff_, var0_, equalvar_, c_] :=
    Block[{var1, var2, stderr},
        {var1, var2} = VariancePair[var0];
        stderr = Sqrt[var1/Length[list1] + var2/Length[list2]];
        NormalCI[diff, stderr, ConfidenceLevel->c]
    ] /; var0 =!= None

VariancePair[{var0_}] := {var0, var0}

VariancePair[var0_List] := {var0[[1]], var0[[2]]} /; Length[var0] >= 2

VariancePair[var0_] := {var0, var0} /; Head[var0] =!= List

(* Ratio of two variances *)

VarianceRatioCI[numlist_List, denlist_List, options___] :=
    Block[{ratio = Variance[numlist]/Variance[denlist],
	   numdof = Length[numlist] - 1, dendof = Length[denlist] - 1,
	   c = (ConfidenceLevel /. {options} /. Options[VarianceRatioCI])},
        FRatioCI[ratio, numdof, dendof, ConfidenceLevel -> c]
    ]

FRatioCI[ratio_, numdof_, dendof_, options___] :=
    Block[{c = ConfidenceLevel /. {options} /. Options[FRatioCI],
	 lower, upper},
        lower = 1/Quantile[FRatioDistribution[numdof, dendof], (1+c)/2];
        upper = 1/Quantile[FRatioDistribution[numdof, dendof], (1-c)/2];
        {lower, upper} ratio
    ]

End[]
SetAttributes[MeanCI ,ReadProtected];
SetAttributes[VarianceCI ,ReadProtected];
SetAttributes[MeanDifferenceCI ,ReadProtected];
SetAttributes[VarianceRatioCI ,ReadProtected];
SetAttributes[StudentTCI ,ReadProtected];
SetAttributes[NormalCI ,ReadProtected];
SetAttributes[ChiSquareCI ,ReadProtected];
SetAttributes[FRatioCI ,ReadProtected];
SetAttributes[ConfidenceLevel ,ReadProtected];
Protect[MeanCI, VarianceCI,MeanDifferenceCI,
        VarianceRatioCI,StudentTCI,
        NormalCI, ChiSquareCI,FRatioCI,ConfidenceLevel];


EndPackage[]


