home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / STATISTI.PAK / HYPOTHES.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  15.4 KB  |  439 lines
  1.  
  2. (*:Version: Mathematica 2.0 *)
  3.  
  4. (*:Name: Statistics`HypothesisTests *)
  5.  
  6. (*:Title: Hypothesis Tests Related to the Normal Distribution *)
  7.  
  8. (*:Author:
  9.   David Withoff (Wolfram Research), February, 1990
  10. *)
  11.  
  12. (*:Legal:
  13.   Copyright (c) 1990, Wolfram Research, Inc.
  14. *)
  15.  
  16. (*:Reference: Usage messages only. *)
  17.  
  18. (*:Summary:
  19.   Hypothesis tests based on elementary distributions derived
  20.   from the normal distribution.  Distributions represented are
  21.   NormalDistribution, StudentTDistribution, ChiSquareDistribution,
  22.   and FRatioDistribution.
  23. *)
  24.  
  25. (*:Keywords: hypothesis test, significance level *)
  26.  
  27. (*:Requirements: No special system requirements. *)
  28.  
  29. (*:Warning: None. *)
  30.  
  31. (*:Sources: Basic statistics texts. *)
  32.  
  33. BeginPackage["Statistics`HypothesisTests`",
  34.          "Statistics`Common`DistributionsCommon`",
  35.              "Statistics`NormalDistribution`",
  36.              "Statistics`DescriptiveStatistics`",
  37.              "Statistics`Common`HypothesisCommon`"]
  38.  
  39. ResultOfTest::usage =
  40. "ResultOfTest[pvalue, options] is a utility used by the hypothesis testing
  41. functions to determine the conclusion of the test using the probability
  42. estimate pvalue for the test described by the options."
  43.  
  44. (* Note: The option list for ResultOfTest must appear in this package   *)
  45. (* before the other test function option lists, since it is used in      *)
  46. (* assigning default values for those lists.                            *)
  47. Options[ResultOfTest] = {SignificanceLevel -> None, TwoSided -> False}
  48.  
  49. MeanTest::usage =
  50. "MeanTest[list, mu0, options] returns a probability estimate (p-value)
  51. for the relationship between the hypothesized population mean mu0 and
  52. the mean of the entries in list."
  53.  
  54. Options[MeanTest] = {KnownStandardDeviation -> None,FullReport -> False,
  55.     KnownVariance -> None} ~Join~ Options[ResultOfTest]
  56.  
  57. MeanDifferenceTest::usage =
  58. "MeanDifferenceTest[list1, list2, diff0, options] returns probability
  59. estimates (p-values) and other hypothesis test information for the
  60. relationship between the hypothesized population mean difference diff0
  61. and Mean[list1] - Mean[list2]."
  62.  
  63. Options[MeanDifferenceTest] = {KnownStandardDeviation -> None,
  64.     FullReport->False, KnownVariance -> None, 
  65.     EqualVariances -> False} ~Join~  Options[ResultOfTest]
  66.  
  67. VarianceTest::usage =
  68. "VarianceTest[list, var0, options] returns probability
  69. estimates (p-values) and other hypothesis test information for
  70. the relationship between var0 and Variance[list]."
  71.  
  72. Options[VarianceTest] = {FullReport->False} ~Join~ Options[ResultOfTest]
  73.  
  74. VarianceRatioTest::usage =
  75. "VarianceRatioTest[numlist, denlist, ratio0, options] returns
  76. probability estimates (p-values) and other hypothesis test information
  77. for the relationship between ratio0 and the ratio
  78. Variance[numlist]/Variance[denlist]."
  79.  
  80. Options[VarianceRatioTest] = {FullReport->False} ~Join~ Options[ResultOfTest]
  81.  
  82. NormalPValue::usage =
  83. "NormalPValue[teststat] returns the cumulative density beyond teststat
  84. for NormalDistribution[]."
  85.  
  86. Options[NormalPValue] = Options[ResultOfTest]
  87.  
  88. StudentTPValue::usage =
  89. "StudentTPValue[teststat, dof] returns the cumulative density beyond
  90. teststat for the StudentTDistribution, with dof degrees of freedom."
  91.  
  92. Options[StudentTPValue] = Options[ResultOfTest]
  93.  
  94. ChiSquarePValue::usage =
  95. "ChiSquarePValue[teststat, dof] returns the cumulative density beyond 
  96. teststat for the ChiSquareDistribution with dof degrees of freedom."
  97.  
  98. Options[ChiSquarePValue] = Options[ResultOfTest]
  99.  
  100. FRatioPValue::usage =
  101. "FRatioPValue[teststat, numdof, dendof] returns the cumulative
  102. density beyond teststat for the FRatioDistribution with numdof numerator 
  103. degrees of freedom and dendof denominator degrees of freedom."
  104.  
  105. Options[FRatioPValue] = Options[ResultOfTest]
  106.  
  107.     (* Output names *)
  108.  
  109. OneSidedPValue::usage =
  110. "OneSidedPValue is used in the output of statistical hypothesis tests
  111. to identify the probability for sample parameters to be further from
  112. the corresponding population parameter than is the current sample parameter,
  113. and on the same side of the sampling distribution."
  114.  
  115. TwoSidedPValue::usage =
  116. "TwoSidedPValue is used in the output of statistical hypothesis tests
  117. to identify the probability for sample parameters to be further from
  118. the corresponding population parameter than the current sample parameter,
  119. on either side of the sampling distribution."
  120.  
  121.     (*  Options *)
  122.  
  123. SignificanceLevel::usage =
  124. "SignificanceLevel is an option to statistical hypothesis tests and is
  125. used to specify the significance level of the test."
  126.  
  127. TwoSided::usage =
  128. "TwoSided is an option to statistical hypothesis tests and is used
  129. to request a two-sided test."
  130.  
  131. FullReport::usage = 
  132. "FullReport is an option to statistical hypothesis tests and is used
  133. to indicate whether the estimator, test statistic, and number of
  134. degrees of freedom are included in the output."
  135.  
  136.  
  137. Begin["`Private`"]
  138.  
  139. VariancePair[{var0_}] := {var0, var0}
  140.  
  141. VariancePair[var0_List] := {var0[[1]], var0[[2]]} /; Length[var0] >= 2
  142.  
  143. VariancePair[var0_] := {var0, var0} /; Head[var0] =!= List
  144.  
  145. (* ==== Hypothesis tests for one mean ================================= *)
  146.  
  147. MeanTest[args___] :=
  148.     Block[{answer = iMeanTest[OptionExtract[{args}, MeanTest]]},
  149.         answer /; answer =!= Fail
  150.     ]
  151.  
  152. iMeanTest[{list_List?(Length[#] > 0 &), mu0_, optionlist_}] :=
  153.     Block[{diff, n, var0, sd0, test, subopts, out, rep},
  154.         diff = Mean[list] - mu0;
  155.         n = Length[list];
  156.         {var0, sd0, rep} = ReplaceAll[{KnownVariance, KnownStandardDeviation,
  157.             FullReport}, optionlist] /. Options[MeanTest];
  158.         If[sd0 =!= None,
  159.             If[var0 =!= None, Message[MeanTest::varsd]];
  160.             var0 = sd0^2
  161.         ];
  162.         If[var0 =!= None,
  163.             test = diff Sqrt[n/var0];
  164.             subopts = OptionSelect[optionlist, NormalPValue];
  165.             out = NormalPValue[test, subopts];
  166.         If[TrueQ[rep],
  167.          rep = FullReport ->TableForm[{{Mean[list],N[test]}},
  168.         TableHeadings->{None,{"Mean","TestStat"}}];
  169.         out =Flatten[{rep,"NormalDistribution",out}]];
  170.         out,    
  171.         (* else, estimate the variance from the sample. *)
  172.             var0 = VarianceOfSampleMean[list];
  173.             If[Head[var0] === VarianceOfSampleMean,
  174.                 Message[MeanTest::novar]; Return[Fail]
  175.             ];
  176.             If[var0 == 0, Message[MeanTest::zerovar]];
  177.             test = diff/Sqrt[var0];
  178.             subopts = OptionSelect[optionlist, StudentTPValue];
  179.             out = StudentTPValue[test, n-1, subopts];
  180.         If[TrueQ[rep],
  181.         rep = FullReport ->TableForm[{{Mean[list],N[test],n-1}},
  182.         TableHeadings->{None,{"Mean","TestStat","DF"}}];
  183.         out =Flatten[{rep,"StudentTDistribution",out}]];
  184.         out   ]]
  185.     
  186.  
  187. iMeanTest[badargs_] :=
  188.     (If[badargs =!= Fail, Message[MeanTest::badargs]]; Fail)
  189.  
  190. MeanTest::badargs = "Incorrect number or type of arguments."
  191.  
  192. MeanTest::varsd = "Warning: KnownStandardDeviation and KnownVariance
  193. have both been specified.  KnownVariance will be ignored."
  194.  
  195. MeanTest::novar = "Unable to estimate variance from the sample."
  196.  
  197. MeanTest::zerovar = "Warning: Estimated variance is zero; subsequent
  198. results may be misleading."
  199.  
  200. (* ==== Hypothesis test for a difference of means ===================== *)
  201.  
  202. MeanDifferenceTest[args___] :=
  203.     Block[{result = vMeanDifferenceTest[
  204.                 OptionExtract[{args}, MeanDifferenceTest]]},
  205.         result /; result =!= Fail]
  206.  
  207. vMeanDifferenceTest[{list1_List, list2_List, diff1minus2_, optionlist_}] :=
  208.     Block[{diff, var0, sd0, pval},
  209.         diff = Mean[list1] - Mean[list2] - diff1minus2;
  210.         {var0, sd0} = ReplaceAll[{KnownVariance, KnownStandardDeviation},
  211.                 optionlist] /. Options[MeanDifferenceTest];
  212.         If[sd0 =!= None,
  213.             If[var0 =!= None, Message[MeanDifferenceTest::varsd]];
  214.             var0 = sd0^2
  215.         ];
  216.         iMeanDifferenceTest[list1, list2, diff, var0, optionlist]
  217.  
  218.     ]
  219.  
  220. vMeanDifferenceTest[badargs_] :=
  221.     (If[badargs =!= Fail, Message[MeanDifferenceTest::badargs]]; Fail)
  222.  
  223. MeanDifferenceTest::badargs = "Incorrect number or type of arguments."
  224.  
  225. MeanDifferenceTest::varsd = "Warning: KnownStandardDeviation and
  226. KnownVariance have both been specified.  KnownVariance will be ignored."
  227.  
  228. iMeanDifferenceTest[list1_, list2_, diff_, None, optionlist_] :=
  229.     Block[
  230.     {var1, var2, n1, n2, dof, pooledvar, test, equalvar, subopts,rep,out},
  231.         {equalvar,rep} = {EqualVariances,FullReport} /. optionlist;
  232.         var1 = Variance[list1];
  233.         var2 = Variance[list2];
  234.         n1 = Length[list1];
  235.         n2 = Length[list2];
  236.         If[TrueQ[equalvar],
  237.             dof = n1 + n2 - 2;
  238.             pooledvar = ((n1-1) var1 + (n2-1) var2) / dof;
  239.             test = diff/Sqrt[pooledvar (1/n1 + 1/n2)],
  240.         (* else *)
  241.             pooledvar = var1/n1 + var2/n2;
  242.             dof = pooledvar^2 /
  243.                       ((var1/n1)^2/(n1-1) + (var2/n2)^2/(n2-1));
  244.             test = diff/Sqrt[pooledvar]
  245.         ];
  246.         subopts = OptionSelect[optionlist, StudentTPValue];
  247.         out = StudentTPValue[test, dof, subopts];
  248.      If[TrueQ[rep],
  249.        rep = FullReport ->TableForm[{{N[Mean[list1]-Mean[list2]],N[test],dof
  250.     }}  ,TableHeadings->{None,{"MeanDiff","TestStat","DF"}}];
  251.     out =Flatten[{rep,"StudentTDistribution",out}]];
  252.         out   
  253.     ]
  254.  
  255. iMeanDifferenceTest[list1_, list2_, diff_, var0_, options___] :=
  256.     Block[{var1, var2, n1, n2, dof, pooledvar, test, subopts,out,rep},
  257.         {var1, var2} = VariancePair[var0];
  258.     rep = FullReport /. options;
  259.         n1 = Length[list1];
  260.         n2 = Length[list2];
  261.         pooledvar = var1/n1 + var2/n2;
  262.         test = diff/Sqrt[pooledvar];
  263.         subopts = OptionSelect[{options}, NormalPValue];
  264.         out = NormalPValue[test, subopts];
  265.     If[TrueQ[rep],
  266.       rep = FullReport ->TableForm[{{N[Mean[list1]-Mean[list2]],N[test]}}
  267.       ,TableHeadings->{None,{"MeanDiff","TestStat"}}];
  268.        out =Flatten[{rep,"NormalDistribution",out}]];
  269.        out   
  270.     ] /; var0 =!= None
  271.  
  272. (* ==== Hypothesis test for a single variance ========================= *)
  273.  
  274. VarianceTest[args___] :=
  275.     Block[{result = iVarianceTest[
  276.                 OptionExtract[{args}, VarianceTest]]},
  277.         result /; result =!= Fail]
  278.  
  279. iVarianceTest[{list_List, var0_, optionlist_}] :=
  280.     Block[{dof, ssq, test, subopts, rep, out},
  281.     rep = FullReport /. optionlist /. Options[VarianceTest];
  282.         dof = Length[list]-1;
  283.         ssq = (dof+1) Variance[list];
  284.         test = ssq/var0;
  285.         subopts = OptionSelect[optionlist, ChiSquarePValue];
  286.         out = ChiSquarePValue[test, dof, subopts];
  287.     If[TrueQ[rep],
  288.         rep = FullReport ->TableForm[{{Variance[list],N[test],dof}},
  289.         TableHeadings->{None,{"Variance","TestStat","DF"}}];
  290.         out =Flatten[{rep,"ChiSquare Distribution",out}]]; out
  291.     ]
  292.  
  293. iVarianceTest[badargs_] :=
  294.     (If[badargs =!= Fail, Message[VarianceTest::badargs]]; Fail)
  295.  
  296. VarianceTest::badargs = "Incorrect number or type of arguments."
  297.  
  298. (* ==== Hypothesis test for a variance ratio ========================== *)
  299.  
  300. VarianceRatioTest[args___] :=
  301.     Block[{result = iVarianceRatioTest[
  302.                 OptionExtract[{args}, VarianceRatioTest]]},
  303.         result /; result =!= Fail]
  304.  
  305. iVarianceRatioTest[
  306.         {numlist_List, denlist_List, ratio0_, optionlist_}] :=
  307.     Block[{numdof, dendof, test, pval, subopts, rep, out},
  308.     rep = FullReport /.optionlist /. Options[VarianceRatioTest];
  309.         test = (Variance[numlist]/Variance[denlist])/ratio0;
  310.         numdof = Length[numlist] - 1;
  311.         dendof = Length[denlist] - 1;
  312.         subopts = OptionSelect[optionlist, FRatioPValue];
  313.         out = FRatioPValue[test, numdof, dendof, subopts];
  314.      If[TrueQ[rep],
  315.         rep=FullReport->TableForm[{{N[test ratio0],N[test],numdof,dendof}}
  316.         ,TableHeadings->{None,{"Ratio","TestStat","NumDF","DenDF"}}];
  317.         out =Flatten[{rep,"FRatio Distribution",out}]]; out
  318.     ]
  319.  
  320. iVarianceRatioTest[badargs_] :=
  321.     CompoundExpression[ 
  322.         If[badargs =!= Fail,
  323.             Message[VarianceRatioTest::badargs, badargs]];
  324.         Fail
  325.     ]
  326.  
  327. VarianceRatioTest::badargs = "Incorrect arguments `1`."
  328.  
  329. (* ==== Basic hypothesis test functions ========================== *)
  330.  
  331. NormalPValue[test_, options___] :=
  332.     Block[{pval},
  333.         pval = N[CDF[NormalDistribution[], test]];
  334.         ResultOfTest[pval, options]
  335.     ]
  336.  
  337. StudentTPValue[test_, dof_, options___] :=
  338.     Block[{pval},
  339.         pval = N[CDF[StudentTDistribution[dof], test]];
  340.         ResultOfTest[pval, options]
  341.     ]
  342.  
  343. ChiSquarePValue[test_, dof_, options___] :=
  344.  Block[{pval},
  345.         pval = N[CDF[ChiSquareDistribution[dof], test]];
  346.         ResultOfTest[pval, options]
  347.     ]
  348.  
  349. FRatioPValue[test_, numdof_, dendof_, options___] :=
  350.     Block[{pval},
  351.         pval = N[CDF[FRatioDistribution[numdof, dendof], test]];
  352.         ResultOfTest[pval, options]
  353.     ] 
  354.  
  355. (*  Hypothesis test utilities  *)
  356.  
  357. ResultOfTest[pval_, options___] :=
  358.     Block[{properpval, sig, twosided, report},
  359.         {sig, twosided} = {SignificanceLevel, TwoSided} /.
  360.             Options[{options}] /. Options[ResultOfTest];
  361.         properpval = ProperPValue[pval];
  362.         If[TrueQ[twosided],
  363.             properpval = 2 properpval;
  364.             report = TwoSidedPValue -> properpval,
  365.         (* else *)
  366.             report = OneSidedPValue -> properpval
  367.         ];
  368.         If[sig =!= None,
  369.             report = {report, SignificanceMessage[properpval, sig]}
  370.         ];
  371.         report
  372.     ]
  373.  
  374. ProperPValue[pval_] :=
  375.     If[pval > .5,
  376.         1 - pval,
  377.     (* else *)
  378.         pval,
  379.     (* undetermined *)
  380.         Message[ResultOfTest::nonum, pval];
  381.         pval
  382.     ]
  383.  
  384. ResultOfTest::nonum = "Warning: P-value `1` is non-numerical.
  385. The symbolic answer may be ambiguous."
  386.  
  387. SignificanceMessage[pval_, level_] :=
  388.     If[N[pval > level],
  389.         "Accept null hypothesis at significance level" -> level,
  390.     (* else *)
  391.         "Reject null hypothesis at significance level" -> level,
  392.     (* undetermined *)
  393.         If[pval > level,
  394.             "Accept null hypothesis at significance level" -> level]
  395.     ]
  396.  
  397. OptionExtract[input_List, f_] :=
  398.     Module[{n, opts, answer, known},
  399.         For[n = Length[input], n > 0, n--,
  400.             If[!OptionQ[input[[n]]], Break[]]
  401.         ];
  402.         answer = Take[input, n];
  403.         opts = Options[input];
  404.         known = Map[First, Options[f]];
  405.         opts = Select[opts,
  406.             If[MemberQ[known,First[#]], True,
  407.                 Message[f::optx, #, f]; False] &];
  408.         AppendTo[answer, opts]
  409.     ]
  410.  
  411. OptionSelect[opts_List, sel_] := 
  412.     With[{known = First /@ Options[sel]}, 
  413.         Select[Flatten[opts], MemberQ[known, First[#]] &]
  414.     ]
  415.  
  416. End[]
  417. SetAttributes[ ResultOfTest ,ReadProtected];
  418. SetAttributes[ MeanTest ,ReadProtected];
  419. SetAttributes[ MeanDifferenceTest ,ReadProtected];
  420. SetAttributes[ VarianceTest ,ReadProtected];
  421. SetAttributes[ VarianceRatioTest ,ReadProtected];
  422. SetAttributes[ NormalPValue ,ReadProtected];
  423. SetAttributes[ StudentTPValue ,ReadProtected];
  424. SetAttributes[ ChiSquarePValue, ReadProtected];
  425. SetAttributes[ FRatioPValue, ReadProtected];
  426. SetAttributes[ OneSidedPValue, ReadProtected];
  427. SetAttributes[ TwoSidedPValue, ReadProtected];
  428. SetAttributes[ SignificanceLevel, ReadProtected];
  429. SetAttributes[ TwoSided, ReadProtected];
  430. SetAttributes[ SignTest, ReadProtected];
  431.  
  432. Protect[ResultOfTest,MeanTest,MeanDifferenceTest,VarianceTest, 
  433.     VarianceRatioTest, NormalPValue,StudentTPValue, ChiSquarePValue,
  434.     FRatioPValue,OneSidedPValue, TwoSidedPValue,SignificanceLevel,
  435.     TwoSided, SignTest];
  436.  
  437.  
  438. EndPackage[]
  439.