(*:Mathematica Version: 2.0 *)

(*:Package Version: 1.1 *)

(*:Name: Statistics`HypothesisTests` *)

(*:Context: Statistics`HypothesisTests` *)

(*:Title: Hypothesis Tests Related to the Normal Distribution *)

(*:Author:
  David Withoff (Wolfram Research), February 1990
*)

(*:History:
  Original version by David Withoff (Wolfram Research), February 1990.
  Changed FullReport->True to give a valid list of replacement rules
	and the distribution in symbolic form;  added LargeSampleSize
	option, and PValueTable symbol;  changed significance message to
	say "fail to reject" rather than "accept"; changed NormalPValue,
	StudentTPValue, ChiSquarePValue, and FRatioPValue to give correct
	results for symbolic input, ECM (Wolfram Research), 1994, 1996.
*)

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

(*:Reference: Usage messages only. *)

(*:Summary:
This package provides hypothesis tests based on elementary distributions
derived from the normal distribution.  Distributions represented are
NormalDistribution, StudentTDistribution, ChiSquareDistribution,
and FRatioDistribution.
*)

(*:Keywords: hypothesis test, significance level *)

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

(*:Warning: None. *)

(*:Sources: Basic statistics texts. *)

BeginPackage["Statistics`HypothesisTests`",
	(* needed for NormalDistribution, StudentTDistribution *)	
             "Statistics`NormalDistribution`",
	(* needed for Mean, Variance *)
             "Statistics`DescriptiveStatistics`",
	(* needed for KnownStandardDeviation, KnownVariance, EqualVariances *)
	     "Statistics`Common`PopulationsCommon`",
	(* needed for CDF *)
	     "Statistics`Common`DistributionsCommon`"]

ResultOfTest::usage =
"ResultOfTest[pvalue, options] gives the conclusion of the test
described by the options using the probability estimate pvalue."

Options[ResultOfTest] = {SignificanceLevel -> None, TwoSided -> False}
 
MeanTest::usage =
"MeanTest[list, mu0, options] returns a probability estimate (p-value)
and other hypthesis test information for the relationship between the
hypothesized population mean mu0 and Mean[list]."

Options[MeanTest] = {KnownStandardDeviation -> None, FullReport -> False,
    KnownVariance -> None} ~Join~ Options[ResultOfTest]

MeanDifferenceTest::usage =
"MeanDifferenceTest[list1, list2, diff0, options] returns a probability
estimate (p-value) and other hypothesis test information for the
relationship between the hypothesized population mean difference diff0
and Mean[list1] - Mean[list2]."

Options[MeanDifferenceTest] = {KnownStandardDeviation -> None,
    FullReport->False, KnownVariance -> None, 
    EqualVariances -> False} ~Join~  Options[ResultOfTest]

VarianceTest::usage =
"VarianceTest[list, var0, options] returns a probability estimate (p-value)
and other hypothesis test information for the relationship between the
hypothesized population variance var0 and Variance[list]."

Options[VarianceTest] = {FullReport->False} ~Join~ Options[ResultOfTest]

VarianceRatioTest::usage =
"VarianceRatioTest[numlist, denlist, ratio0, options] returns a probability
estimate (p-value) and other hypothesis test information for the relationship
between the hypothesized population variance ratio ratio0 and the the ratio
Variance[numlist]/Variance[denlist]."

Options[VarianceRatioTest] = {FullReport->False} ~Join~ Options[ResultOfTest]

NormalPValue::usage =
"NormalPValue[teststat] returns the cumulative density beyond teststat
for NormalDistribution[]."

Options[NormalPValue] = Options[ResultOfTest]

StudentTPValue::usage =
"StudentTPValue[teststat, dof] returns the cumulative density beyond
teststat for the StudentTDistribution, with dof degrees of freedom."

Options[StudentTPValue] = Options[ResultOfTest]

ChiSquarePValue::usage =
"ChiSquarePValue[teststat, dof] returns the cumulative density beyond 
teststat for the ChiSquareDistribution with dof degrees of freedom."

Options[ChiSquarePValue] = Options[ResultOfTest]

FRatioPValue::usage =
"FRatioPValue[teststat, numdof, dendof] returns the cumulative
density beyond teststat for the FRatioDistribution with numdof numerator 
degrees of freedom and dendof denominator degrees of freedom."

Options[FRatioPValue] = Options[ResultOfTest]

	(* Output names *)

OneSidedPValue::usage =
"OneSidedPValue is used in the output of statistical hypothesis tests to
identify the probability of observing a value further from the population
parameter than is the sample parameter, and on the same side of the sampling
distribution."

TwoSidedPValue::usage =
"TwoSidedPValue is used in the output of statistical hypothesis tests to
identify the probability of observing a value further from the population
parameter than is the sample parameter, on either side of the sampling
distribution."

PValueTable::usage =
"PValueTable is used in the output of statistical hypothesis tests to
identify a table of p-values associated with a multiple comparison."

	(*  Options *)

SignificanceLevel::usage =
"SignificanceLevel is an option to statistical hypothesis tests and is
used to specify the significance level of the test."

TwoSided::usage =
"TwoSided is an option to statistical hypothesis tests and is used
to request a two-sided test."

FullReport::usage = 
"FullReport is an option to statistical hypothesis tests and is used
to indicate whether such information as the estimator, test statistic, and
number of degrees of freedom should be included in the output."

LargeSampleSize::usage =
"LargeSampleSize is an option for hypothesis testing functions that return
exact p-values.  LargeSampleSize -> n specifies that for a data sample larger
than n, a large-sample approximation to the p-value will be given."


Begin["`Private`"]

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

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

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

(* ==== Hypothesis tests for one mean ================================= *)

MeanTest[args___] :=
    Block[{answer = iMeanTest[OptionExtract[{args}, MeanTest]]},
        answer /; answer =!= Fail
    ]

iMeanTest[{list_List?(Length[#] > 0 &), mu0_, optionlist_}] :=
    Block[{mean = Mean[list], delta, n, var0, sd0, rep, test, subopts,
		 out, nmean, ntest},
        delta = mean - mu0;
        n = Length[list];
        {var0, sd0, rep} = ReplaceAll[{KnownVariance, KnownStandardDeviation,
            FullReport}, optionlist] /. Options[MeanTest];
        If[sd0 =!= None,
            If[var0 =!= None, Message[MeanTest::varsd]];
            var0 = sd0^2
        ];
	nmean = If[Precision[mean] === Infinity, N[mean], mean];
        If[var0 =!= None,
            test = delta/Sqrt[var0/n];
	    ntest = If[Precision[test] === Infinity, N[test], test];  
            subopts = OptionSelect[optionlist, NormalPValue];
            out = NormalPValue[ntest, subopts];
	    If[TrueQ[rep],
	     (* NOTE: 8/94: changed full report to give a valid list of
		 replacement rules and the distribution in symbolic form *)	
	     rep = FullReport ->
		 TableForm[{{nmean, ntest, NormalDistribution[]}},
		TableHeadings -> {None, {"Mean","TestStat", "Distribution"}}];
	     out = Flatten[{rep, out}]
	    ];
	    out,	
        (* else, estimate the variance from the sample. *)
            var0 = VarianceOfSampleMean[list];
            If[Head[var0] === VarianceOfSampleMean,
                Message[MeanTest::novar]; Return[Fail] ];
            If[var0 == 0, Message[MeanTest::zerovar]];
            test = delta/Sqrt[var0];
	    ntest = If[Precision[test] === Infinity, N[test], test];
            subopts = OptionSelect[optionlist, StudentTPValue];
            out = StudentTPValue[ntest, n-1, subopts];
	    If[TrueQ[rep],
	        (* NOTE: 8/94: changed full report to give a valid list of
		   replacement rules and the distribution in symbolic form *)	
		rep = FullReport ->
		 TableForm[{{nmean, ntest, StudentTDistribution[n-1]}},
		 TableHeadings -> {None, {"Mean", "TestStat", "Distribution"}}];
		out = Flatten[{rep, out}]
	    ];
	    out
	]
    ]
    

iMeanTest[badargs_] :=
    (If[badargs =!= Fail, Message[MeanTest::badargs]]; Fail)

MeanTest::badargs = "Incorrect number or type of arguments."

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

MeanTest::novar = "Unable to estimate variance from the sample."

MeanTest::zerovar = "Warning: Estimated variance is zero; subsequent
results may be misleading."

(* ==== Hypothesis test for a difference of means ===================== *)

MeanDifferenceTest[args___] :=
    Block[{result = vMeanDifferenceTest[
                OptionExtract[{args}, MeanDifferenceTest]]},
        result /; result =!= Fail]

vMeanDifferenceTest[{list1_List, list2_List, diff1minus2_, optionlist_}] :=
    Block[{meandiff = Mean[list1] - Mean[list2], diff, var0, sd0, pval},
        delta = meandiff - diff1minus2;
        {var0, sd0} = ReplaceAll[{KnownVariance, KnownStandardDeviation},
                optionlist] /. Options[MeanDifferenceTest];
        If[sd0 =!= None,
            If[var0 =!= None, Message[MeanDifferenceTest::varsd]];
            var0 = sd0^2
        ];
        iMeanDifferenceTest[list1, list2, delta, var0, optionlist]

    ]

vMeanDifferenceTest[badargs_] :=
    (If[badargs =!= Fail, Message[MeanDifferenceTest::badargs]]; Fail)

MeanDifferenceTest::badargs = "Incorrect number or type of arguments."

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

iMeanDifferenceTest[list1_, list2_, delta_, None, optionlist_] :=
    Block[{equalvar, rep, var1, var2, n1, n2, dof, pooledvar, test, subopts,
		 out, meandiff = Mean[list1]-Mean[list2], nmeandiff, ntest},
        {equalvar, rep} = {EqualVariances, FullReport} /. optionlist;
        {var1, var2} = {Variance[list1], Variance[list2]};
        {n1, n2} = {Length[list1], Length[list2]};
        If[TrueQ[equalvar],
            dof = n1 + n2 - 2;
            pooledvar = ((n1-1) var1 + (n2-1) var2) / dof;
            test = delta/Sqrt[pooledvar (1/n1 + 1/n2)],
        (* else *)
            pooledvar = var1/n1 + var2/n2;
            dof = pooledvar^2 /
                      ((var1/n1)^2/(n1-1) + (var2/n2)^2/(n2-1));
            test = delta/Sqrt[pooledvar]
        ];
	ntest = If[Precision[test] === Infinity, N[test], test];
        subopts = OptionSelect[optionlist, StudentTPValue];
        out = StudentTPValue[ntest, dof, subopts];
        If[TrueQ[rep],
	   (* NOTE: 8/94: changed full report to give a valid list of
		 replacement rules and the distribution in symbolic form *)	
	   nmeandiff =
                 If[Precision[meandiff] === Infinity, N[meandiff], meandiff];
           rep = FullReport ->
		 TableForm[{{nmeandiff, ntest, StudentTDistribution[dof]}},
		 TableHeadings -> {None,
			{"MeanDiff", "TestStat", "Distribution"}}];
	   out = Flatten[{rep, out}]
	];
	out   
    ]

iMeanDifferenceTest[list1_, list2_, delta_, var0_, options___] :=
    Block[{rep, var1, var2, n1, n2, pooledvar, test, subopts, out,
		meandiff = Mean[list1]-Mean[list2], nmeandiff, ntest},
	rep = FullReport /. options;
        {var1, var2} = VariancePair[var0];
        {n1, n2} = {Length[list1], Length[list2]};
        pooledvar = var1/n1 + var2/n2;
        test = delta/Sqrt[pooledvar];
	ntest = If[Precision[test] === Infinity, N[test], test]; 
        subopts = OptionSelect[{options}, NormalPValue];
        out = NormalPValue[ntest, subopts];
	If[TrueQ[rep],
	  (* NOTE: 8/94: changed full report to give a valid list of
		 replacement rules and the distribution in symbolic form *)	
	  nmeandiff =
		 If[Precision[meandiff] === Infinity, N[meandiff], meandiff];
	  rep = FullReport ->
		 TableForm[{{nmeandiff, ntest, NormalDistribution[]}},
	  	 TableHeadings -> {None,
			{"MeanDiff","TestStat", "Distribution"}}];
	  out = Flatten[{rep, out}]
        ];
	out   
    ] /; var0 =!= None

(* ==== Hypothesis test for a single variance ========================= *)

VarianceTest[args___] :=
    Block[{result = iVarianceTest[
                OptionExtract[{args}, VarianceTest]]},
        result /; result =!= Fail]

iVarianceTest[{list_List, var0_, optionlist_}] :=
  Block[{rep, dof, ssq, test, subopts, out, variance = Variance[list],
		nvariance, ntest},
    rep = FullReport /. optionlist /. Options[VarianceTest];
    dof = Length[list]-1;
    ssq = dof variance; 
    test = ssq/var0;
    ntest = If[Precision[test] === Infinity, N[test], test];
    subopts = OptionSelect[optionlist, ChiSquarePValue];
    out = ChiSquarePValue[ntest, dof, subopts];
    If[TrueQ[rep],
       (* NOTE: 8/94: changed full report to give a valid list of
		 replacement rules and the distribution in symbolic form *)	
       nvariance = If[Precision[variance] === Infinity, N[variance], variance];	
       rep = FullReport ->
	  TableForm[{{nvariance, ntest, ChiSquareDistribution[dof]}},
	   TableHeadings -> {None, {"Variance", "TestStat", "Distribution"}}];
       out = Flatten[{rep, out}]	];
    out
  ]

iVarianceTest[badargs_] :=
    (If[badargs =!= Fail, Message[VarianceTest::badargs]]; Fail)

VarianceTest::badargs = "Incorrect number or type of arguments."

(* ==== Hypothesis test for a variance ratio ========================== *)

VarianceRatioTest[args___] :=
    Block[{result = iVarianceRatioTest[
                OptionExtract[{args}, VarianceRatioTest]]},
        result /; result =!= Fail]

iVarianceRatioTest[
        {numlist_List, denlist_List, ratio0_, optionlist_}] :=
    Block[{rep, test, numdof, dendof, subopts, out,
	   varRat = Variance[numlist]/Variance[denlist], nvarRat, ntest},
	rep = FullReport /.optionlist /. Options[VarianceRatioTest];
        test = varRat/ratio0;
	ntest = If[Precision[test] === Infinity, N[test], test];
        numdof = Length[numlist] - 1;
        dendof = Length[denlist] - 1;
        subopts = OptionSelect[optionlist, FRatioPValue];
        out = FRatioPValue[ntest, numdof, dendof, subopts];
     If[TrueQ[rep],
	(* NOTE: 8/94: changed full report to give a valid list of
		 replacement rules and the distribution in symbolic form *)	
	nvarRat = If[Precision[varRat] === Infinity, N[varRat], varRat];
        rep = FullReport ->
	  TableForm[{{nvarRat, ntest, FRatioDistribution[numdof, dendof]}},
          TableHeadings -> {None, {"Ratio", "TestStat", "Distribution"}}];
        out = Flatten[{rep, out}]
     ];
     out
    ]

iVarianceRatioTest[badargs_] :=
    CompoundExpression[ 
        If[badargs =!= Fail,
            Message[VarianceRatioTest::badargs, badargs]];
        Fail
    ]

VarianceRatioTest::badargs = "Incorrect arguments `1`."

(* ==== Basic hypothesis test functions ========================== *)

NormalPValue[test_, options___] :=
    Block[{pval = CDF[NormalDistribution[], -Abs[test] ],
	   sig, twosided},
	If[Precision[pval] === Infinity && NumberQ[N[pval]],
	   pval = N[pval]];
        {sig, twosided} = {SignificanceLevel, TwoSided} /.
            Options[{options}] /. Options[NormalPValue];
	resultOfTest[pval, 
                SignificanceLevel -> sig, TwoSided -> twosided]	
    ]

StudentTPValue[test_, dof_, options___] :=
    Block[{pval = CDF[StudentTDistribution[dof], -Abs[test] ],
	   sig, twosided},
	If[Precision[pval] === Infinity && NumberQ[N[pval]],
	   pval = N[pval]];
	(* assume that a symbolic test value is real *)
	pval = pval /. {Sign[Abs[_]] -> 1, Power[Abs[x_], 2] -> x^2};
	{sig, twosided} = {SignificanceLevel, TwoSided} /.
            Options[{options}] /. Options[StudentTPValue];
	resultOfTest[pval, 
                SignificanceLevel -> sig, TwoSided -> twosided]
    ]

ChiSquarePValue[test_, dof_, options___] :=
    Block[{pval = 1/2 - Abs[1/2 - CDF[ChiSquareDistribution[dof], test]],
	   sig, twosided},
	If[Precision[pval] === Infinity && NumberQ[N[pval]],
	   pval = N[pval]];
	{sig, twosided} = {SignificanceLevel, TwoSided} /.
	   Options[{options}] /. Options[ChiSquarePValue];
	resultOfTest[pval, 
                SignificanceLevel -> sig, TwoSided -> twosided]
    ]

FRatioPValue[test_, numdof_, dendof_, options___] :=
    Block[{pval = 1/2 -
		 Abs[1/2 - CDF[FRatioDistribution[numdof, dendof], test]],
	   sig, twosided},
	If[Precision[pval] === Infinity && NumberQ[N[pval]],
	   pval = N[pval]];
        {sig, twosided} = {SignificanceLevel, TwoSided} /.
           Options[{options}] /. Options[FRatioPValue];
	resultOfTest[pval, 
		SignificanceLevel -> sig, TwoSided -> twosided]
    ] 

(*  Hypothesis test utilities  *)

resultOfTest[pval_, options___] :=
    Block[{sig, twosided, report, pval1 = pval},
        {sig, twosided} = {SignificanceLevel, TwoSided} /. Options[{options}];
        If[TrueQ[twosided],
		pval1 = 2 pval1;
		report = TwoSidedPValue -> pval1,
		report = OneSidedPValue -> pval1
	];
        If[sig =!= None,
            report = {report, SignificanceMessage[pval1, sig]}
        ];
        report
    ]

ResultOfTest[pval_, options___] :=
    Block[{properpval = ProperPValue[pval], sig, twosided, report},
        {sig, twosided} = {SignificanceLevel, TwoSided} /.
            Options[{options}] /. Options[ResultOfTest];
	resultOfTest[properpval,
		 SignificanceLevel -> sig, TwoSided -> twosided]
    ]

ProperPValue[pval_] :=
    If[pval > .5,
        1 - pval,
        pval,
        (* undetermined *)
        Message[ResultOfTest::nonum, pval];
        pval
    ]

ResultOfTest::nonum = "Warning: P-value `1` is non-numerical.
The symbolic answer may be ambiguous."
 
SignificanceMessage[pval_, level_] :=
    If[N[pval > level],
        "Fail to reject null hypothesis at significance level" -> level,
        "Reject null hypothesis at significance level" -> level,
         If[pval > level,
        	"Fail to reject null hypothesis at significance level" -> level,
        	"Reject null hypothesis at significance level" -> level
	 ]
    ]

OptionExtract[input_List, f_] :=
    Module[{n, opts, answer, known},
        For[n = Length[input], n > 0, n--,
            If[!OptionQ[input[[n]]], Break[]]
        ];
        answer = Take[input, n];
        opts = Options[input];
        known = Map[First, Options[f]];
        opts = Select[opts,
            If[MemberQ[known,First[#]], True,
                Message[f::optx, #, f]; False] &];
        AppendTo[answer, opts]
    ]

OptionSelect[opts_List, sel_] := 
    With[{known = First /@ Options[sel]}, 
        Select[Flatten[opts], MemberQ[known, First[#]] &]
    ]
 
End[]
SetAttributes[ ResultOfTest ,ReadProtected];
SetAttributes[ MeanTest ,ReadProtected];
SetAttributes[ MeanDifferenceTest ,ReadProtected];
SetAttributes[ VarianceTest ,ReadProtected];
SetAttributes[ VarianceRatioTest ,ReadProtected];
SetAttributes[ NormalPValue ,ReadProtected];
SetAttributes[ StudentTPValue ,ReadProtected];
SetAttributes[ ChiSquarePValue, ReadProtected];
SetAttributes[ FRatioPValue, ReadProtected];
SetAttributes[ OneSidedPValue, ReadProtected];
SetAttributes[ TwoSidedPValue, ReadProtected];
SetAttributes[ SignificanceLevel, ReadProtected];
SetAttributes[ TwoSided, ReadProtected];

Protect[ResultOfTest,MeanTest,MeanDifferenceTest,VarianceTest, 
	VarianceRatioTest, NormalPValue,StudentTPValue, ChiSquarePValue,
	FRatioPValue,OneSidedPValue, TwoSidedPValue,SignificanceLevel,
	TwoSided];


EndPackage[]
