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