f_dist_example.qbk 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220
  1. [section:f_eg F Distribution Examples]
  2. Imagine that you want to compare the standard deviations of two
  3. sample to determine if they differ in any significant way, in this
  4. situation you use the F distribution and perform an F-test. This
  5. situation commonly occurs when conducting a process change comparison:
  6. "is a new process more consistent that the old one?".
  7. In this example we'll be using the data for ceramic strength from
  8. [@http://www.itl.nist.gov/div898/handbook/eda/section4/eda42a1.htm
  9. http://www.itl.nist.gov/div898/handbook/eda/section4/eda42a1.htm].
  10. The data for this case study were collected by Said Jahanmir of the
  11. NIST Ceramics Division in 1996 in connection with a NIST/industry
  12. ceramics consortium for strength optimization of ceramic strength.
  13. The example program is [@../../example/f_test.cpp f_test.cpp],
  14. program output has been deliberately made as similar as possible
  15. to the DATAPLOT output in the corresponding
  16. [@http://www.itl.nist.gov/div898/handbook/eda/section3/eda359.htm
  17. NIST EngineeringStatistics Handbook example].
  18. We'll begin by defining the procedure to conduct the test:
  19. void f_test(
  20. double sd1, // Sample 1 std deviation
  21. double sd2, // Sample 2 std deviation
  22. double N1, // Sample 1 size
  23. double N2, // Sample 2 size
  24. double alpha) // Significance level
  25. {
  26. The procedure begins by printing out a summary of our input data:
  27. using namespace std;
  28. using namespace boost::math;
  29. // Print header:
  30. cout <<
  31. "____________________________________\n"
  32. "F test for equal standard deviations\n"
  33. "____________________________________\n\n";
  34. cout << setprecision(5);
  35. cout << "Sample 1:\n";
  36. cout << setw(55) << left << "Number of Observations" << "= " << N1 << "\n";
  37. cout << setw(55) << left << "Sample Standard Deviation" << "= " << sd1 << "\n\n";
  38. cout << "Sample 2:\n";
  39. cout << setw(55) << left << "Number of Observations" << "= " << N2 << "\n";
  40. cout << setw(55) << left << "Sample Standard Deviation" << "= " << sd2 << "\n\n";
  41. The test statistic for an F-test is simply the ratio of the square of
  42. the two standard deviations:
  43. [expression F = s[sub 1][super 2] / s[sub 2][super 2]]
  44. where s[sub 1] is the standard deviation of the first sample and s[sub 2]
  45. is the standard deviation of the second sample. Or in code:
  46. double F = (sd1 / sd2);
  47. F *= F;
  48. cout << setw(55) << left << "Test Statistic" << "= " << F << "\n\n";
  49. At this point a word of caution: the F distribution is asymmetric,
  50. so we have to be careful how we compute the tests, the following table
  51. summarises the options available:
  52. [table
  53. [[Hypothesis][Test]]
  54. [[The null-hypothesis: there is no difference in standard deviations (two sided test)]
  55. [Reject if F <= F[sub (1-alpha/2; N1-1, N2-1)] or F >= F[sub (alpha/2; N1-1, N2-1)] ]]
  56. [[The alternative hypothesis: there is a difference in means (two sided test)]
  57. [Reject if F[sub (1-alpha/2; N1-1, N2-1)] <= F <= F[sub (alpha/2; N1-1, N2-1)] ]]
  58. [[The alternative hypothesis: Standard deviation of sample 1 is greater
  59. than that of sample 2]
  60. [Reject if F < F[sub (alpha; N1-1, N2-1)] ]]
  61. [[The alternative hypothesis: Standard deviation of sample 1 is less
  62. than that of sample 2]
  63. [Reject if F > F[sub (1-alpha; N1-1, N2-1)] ]]
  64. ]
  65. Where F[sub (1-alpha; N1-1, N2-1)] is the lower critical value of the F distribution
  66. with degrees of freedom N1-1 and N2-1, and F[sub (alpha; N1-1, N2-1)] is the upper
  67. critical value of the F distribution with degrees of freedom N1-1 and N2-1.
  68. The upper and lower critical values can be computed using the quantile function:
  69. [expression F[sub (1-alpha; N1-1, N2-1)] = `quantile(fisher_f(N1-1, N2-1), alpha)`]
  70. [expression F[sub (alpha; N1-1, N2-1)] = `quantile(complement(fisher_f(N1-1, N2-1), alpha))`]
  71. In our example program we need both upper and lower critical values for alpha
  72. and for alpha/2:
  73. double ucv = quantile(complement(dist, alpha));
  74. double ucv2 = quantile(complement(dist, alpha / 2));
  75. double lcv = quantile(dist, alpha);
  76. double lcv2 = quantile(dist, alpha / 2);
  77. cout << setw(55) << left << "Upper Critical Value at alpha: " << "= "
  78. << setprecision(3) << scientific << ucv << "\n";
  79. cout << setw(55) << left << "Upper Critical Value at alpha/2: " << "= "
  80. << setprecision(3) << scientific << ucv2 << "\n";
  81. cout << setw(55) << left << "Lower Critical Value at alpha: " << "= "
  82. << setprecision(3) << scientific << lcv << "\n";
  83. cout << setw(55) << left << "Lower Critical Value at alpha/2: " << "= "
  84. << setprecision(3) << scientific << lcv2 << "\n\n";
  85. The final step is to perform the comparisons given above, and print
  86. out whether the hypothesis is rejected or not:
  87. cout << setw(55) << left <<
  88. "Results for Alternative Hypothesis and alpha" << "= "
  89. << setprecision(4) << fixed << alpha << "\n\n";
  90. cout << "Alternative Hypothesis Conclusion\n";
  91. cout << "Standard deviations are unequal (two sided test) ";
  92. if((ucv2 < F) || (lcv2 > F))
  93. cout << "ACCEPTED\n";
  94. else
  95. cout << "REJECTED\n";
  96. cout << "Standard deviation 1 is less than standard deviation 2 ";
  97. if(lcv > F)
  98. cout << "ACCEPTED\n";
  99. else
  100. cout << "REJECTED\n";
  101. cout << "Standard deviation 1 is greater than standard deviation 2 ";
  102. if(ucv < F)
  103. cout << "ACCEPTED\n";
  104. else
  105. cout << "REJECTED\n";
  106. cout << endl << endl;
  107. Using the ceramic strength data as an example we get the following
  108. output:
  109. [pre
  110. '''F test for equal standard deviations
  111. ____________________________________
  112. Sample 1:
  113. Number of Observations = 240
  114. Sample Standard Deviation = 65.549
  115. Sample 2:
  116. Number of Observations = 240
  117. Sample Standard Deviation = 61.854
  118. Test Statistic = 1.123
  119. CDF of test statistic: = 8.148e-001
  120. Upper Critical Value at alpha: = 1.238e+000
  121. Upper Critical Value at alpha/2: = 1.289e+000
  122. Lower Critical Value at alpha: = 8.080e-001
  123. Lower Critical Value at alpha/2: = 7.756e-001
  124. Results for Alternative Hypothesis and alpha = 0.0500
  125. Alternative Hypothesis Conclusion
  126. Standard deviations are unequal (two sided test) REJECTED
  127. Standard deviation 1 is less than standard deviation 2 REJECTED
  128. Standard deviation 1 is greater than standard deviation 2 REJECTED'''
  129. ]
  130. In this case we are unable to reject the null-hypothesis, and must instead
  131. reject the alternative hypothesis.
  132. By contrast let's see what happens when we use some different
  133. [@http://www.itl.nist.gov/div898/handbook/prc/section3/prc32.htm
  134. sample data]:, once again from the NIST Engineering Statistics Handbook:
  135. A new procedure to assemble a device is introduced and tested for
  136. possible improvement in time of assembly. The question being addressed
  137. is whether the standard deviation of the new assembly process (sample 2) is
  138. better (i.e., smaller) than the standard deviation for the old assembly
  139. process (sample 1).
  140. [pre
  141. '''____________________________________
  142. F test for equal standard deviations
  143. ____________________________________
  144. Sample 1:
  145. Number of Observations = 11.00000
  146. Sample Standard Deviation = 4.90820
  147. Sample 2:
  148. Number of Observations = 9.00000
  149. Sample Standard Deviation = 2.58740
  150. Test Statistic = 3.59847
  151. CDF of test statistic: = 9.589e-001
  152. Upper Critical Value at alpha: = 3.347e+000
  153. Upper Critical Value at alpha/2: = 4.295e+000
  154. Lower Critical Value at alpha: = 3.256e-001
  155. Lower Critical Value at alpha/2: = 2.594e-001
  156. Results for Alternative Hypothesis and alpha = 0.0500
  157. Alternative Hypothesis Conclusion
  158. Standard deviations are unequal (two sided test) REJECTED
  159. Standard deviation 1 is less than standard deviation 2 REJECTED
  160. Standard deviation 1 is greater than standard deviation 2 ACCEPTED'''
  161. ]
  162. In this case we take our null hypothesis as "standard deviation 1 is
  163. less than or equal to standard deviation 2", since this represents the "no change"
  164. situation. So we want to compare the upper critical value at /alpha/
  165. (a one sided test) with the test statistic, and since 3.35 < 3.6 this
  166. hypothesis must be rejected. We therefore conclude that there is a change
  167. for the better in our standard deviation.
  168. [endsect][/section:f_eg F Distribution]
  169. [/
  170. Copyright 2006 John Maddock and Paul A. Bristow.
  171. Distributed under the Boost Software License, Version 1.0.
  172. (See accompanying file LICENSE_1_0.txt or copy at
  173. http://www.boost.org/LICENSE_1_0.txt).
  174. ]