
(*:Version: Mathematica 2.0 *)

(*:Name: Statistics`ConfidenceIntervals *)

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

(*:Author:
  Dave Withoff (Wolfram Research), February, 1990.
  Revised by March, 1991.
*)

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

(*:Reference: usage messages only. *)

(*:Summary:
  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`HypothesisCommon`"]
 
MeanCI::usage =
"MeanCI[list, options] returns a list {min, max}
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 {min, max} 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 {min, max}
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 {min, max} 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
{min, max} 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
{min, max} 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 {min, max} 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, n, c, var0, sd0},
        mean = Mean[list];
        n = Length[list];
        {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, options]],
            Return[NormalCI[mean, Sqrt[var0/n], options]]
        ]
    ]

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, var,
         c = ConfidenceLevel /. {options} /.
             Options[VarianceCI]
        },
        dof = Length[list] - 1;
        var = Variance[list];
        ChiSquareCI[var, dof, options]
    ]

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

(* Comparison of two Means *)

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

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

mdInterval[list1_, list2_, diff_, None, options___] :=
    Block[
        {var1, var2, dof, stderr, n1, n2,
         equalvar = EqualVariances /. {options} /.
             Options[MeanDifferenceCI]
        },
        var1 = Variance[list1];
        var2 = Variance[list2];
        n1 = Length[list1];
        n2 = Length[list2];
        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, options]
    ]

mdInterval[list1_, list2_, diff_, var0_, options___] :=
    Block[{var1, var2, stderr},
        {var1, var2} = VariancePair[var0];
        stderr = Sqrt[var1/Length[list1] + var2/Length[list2]];
        NormalCI[diff, stderr, options]
    ] /; 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, numdof, dendof,
         c = ConfidenceLevel /. {options} /.
             Options[VarianceRatioCI]
        },
        ratio = Variance[numlist]/Variance[denlist];
        numdof = Length[numlist] - 1;
        dendof = Length[denlist] - 1;
        FRatioCI[ratio, numdof, dendof, options]
    ]

FRatioCI[ratio_, numdof_, dendof_, options___] :=
    Block[
        {c = ConfidenceLevel /. {options} /. Options[FRatioCI],
         min, max
        },
        min = 1/Quantile[FRatioDistribution[numdof, dendof], (1+c)/2];
        max = 1/Quantile[FRatioDistribution[numdof, dendof], (1-c)/2];
        {min, max} 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[]

