1 // Copyright John Maddock 2006, 2007
2 // Copyright Paul A. Bristow 2010
3
4 // Use, modification and distribution are subject to the
5 // Boost Software License, Version 1.0.
6 // (See accompanying file LICENSE_1_0.txt
7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
8
9 #include <iostream>
10 using std::cout; using std::endl;
11 using std::left; using std::fixed; using std::right; using std::scientific;
12 #include <iomanip>
13 using std::setw;
14 using std::setprecision;
15
16 #include <boost/math/distributions/chi_squared.hpp>
17
18 int error_result = 0;
19
confidence_limits_on_std_deviation(double Sd,unsigned N)20 void confidence_limits_on_std_deviation(
21 double Sd, // Sample Standard Deviation
22 unsigned N) // Sample size
23 {
24 // Calculate confidence intervals for the standard deviation.
25 // For example if we set the confidence limit to
26 // 0.95, we know that if we repeat the sampling
27 // 100 times, then we expect that the true standard deviation
28 // will be between out limits on 95 occations.
29 // Note: this is not the same as saying a 95%
30 // confidence interval means that there is a 95%
31 // probability that the interval contains the true standard deviation.
32 // The interval computed from a given sample either
33 // contains the true standard deviation or it does not.
34 // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda358.htm
35
36 // using namespace boost::math; // potential name ambiguity with std <random>
37 using boost::math::chi_squared;
38 using boost::math::quantile;
39 using boost::math::complement;
40
41 // Print out general info:
42 cout <<
43 "________________________________________________\n"
44 "2-Sided Confidence Limits For Standard Deviation\n"
45 "________________________________________________\n\n";
46 cout << setprecision(7);
47 cout << setw(40) << left << "Number of Observations" << "= " << N << "\n";
48 cout << setw(40) << left << "Standard Deviation" << "= " << Sd << "\n";
49 //
50 // Define a table of significance/risk levels:
51 double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
52 //
53 // Start by declaring the distribution we'll need:
54 chi_squared dist(N - 1);
55 //
56 // Print table header:
57 //
58 cout << "\n\n"
59 "_____________________________________________\n"
60 "Confidence Lower Upper\n"
61 " Value (%) Limit Limit\n"
62 "_____________________________________________\n";
63 //
64 // Now print out the data for the table rows.
65 for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
66 {
67 // Confidence value:
68 cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]);
69 // Calculate limits:
70 double lower_limit = sqrt((N - 1) * Sd * Sd / quantile(complement(dist, alpha[i] / 2)));
71 double upper_limit = sqrt((N - 1) * Sd * Sd / quantile(dist, alpha[i] / 2));
72 // Print Limits:
73 cout << fixed << setprecision(5) << setw(15) << right << lower_limit;
74 cout << fixed << setprecision(5) << setw(15) << right << upper_limit << endl;
75 }
76 cout << endl;
77 } // void confidence_limits_on_std_deviation
78
confidence_limits_on_std_deviation_alpha(double Sd,double alpha)79 void confidence_limits_on_std_deviation_alpha(
80 double Sd, // Sample Standard Deviation
81 double alpha // confidence
82 )
83 { // Calculate confidence intervals for the standard deviation.
84 // for the alpha parameter, for a range number of observations,
85 // from a mere 2 up to a million.
86 // O. L. Davies, Statistical Methods in Research and Production, ISBN 0 05 002437 X,
87 // 4.33 Page 68, Table H, pp 452 459.
88
89 // using namespace std;
90 // using namespace boost::math;
91 using boost::math::chi_squared;
92 using boost::math::quantile;
93 using boost::math::complement;
94
95 // Define a table of numbers of observations:
96 unsigned int obs[] = {2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30, 40 , 50, 60, 100, 120, 1000, 10000, 50000, 100000, 1000000};
97
98 cout << // Print out heading:
99 "________________________________________________\n"
100 "2-Sided Confidence Limits For Standard Deviation\n"
101 "________________________________________________\n\n";
102 cout << setprecision(7);
103 cout << setw(40) << left << "Confidence level (two-sided) " << "= " << alpha << "\n";
104 cout << setw(40) << left << "Standard Deviation" << "= " << Sd << "\n";
105
106 cout << "\n\n" // Print table header:
107 "_____________________________________________\n"
108 "Observations Lower Upper\n"
109 " Limit Limit\n"
110 "_____________________________________________\n";
111 for(unsigned i = 0; i < sizeof(obs)/sizeof(obs[0]); ++i)
112 {
113 unsigned int N = obs[i]; // Observations
114 // Start by declaring the distribution with the appropriate :
115 chi_squared dist(N - 1);
116
117 // Now print out the data for the table row.
118 cout << fixed << setprecision(3) << setw(10) << right << N;
119 // Calculate limits: (alpha /2 because it is a two-sided (upper and lower limit) test.
120 double lower_limit = sqrt((N - 1) * Sd * Sd / quantile(complement(dist, alpha / 2)));
121 double upper_limit = sqrt((N - 1) * Sd * Sd / quantile(dist, alpha / 2));
122 // Print Limits:
123 cout << fixed << setprecision(4) << setw(15) << right << lower_limit;
124 cout << fixed << setprecision(4) << setw(15) << right << upper_limit << endl;
125 }
126 cout << endl;
127 }// void confidence_limits_on_std_deviation_alpha
128
chi_squared_test(double Sd,double D,unsigned N,double alpha)129 void chi_squared_test(
130 double Sd, // Sample std deviation
131 double D, // True std deviation
132 unsigned N, // Sample size
133 double alpha) // Significance level
134 {
135 //
136 // A Chi Squared test applied to a single set of data.
137 // We are testing the null hypothesis that the true
138 // standard deviation of the sample is D, and that any variation is down
139 // to chance. We can also test the alternative hypothesis
140 // that any difference is not down to chance.
141 // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda358.htm
142 //
143 // using namespace boost::math;
144 using boost::math::chi_squared;
145 using boost::math::quantile;
146 using boost::math::complement;
147 using boost::math::cdf;
148
149 // Print header:
150 cout <<
151 "______________________________________________\n"
152 "Chi Squared test for sample standard deviation\n"
153 "______________________________________________\n\n";
154 cout << setprecision(5);
155 cout << setw(55) << left << "Number of Observations" << "= " << N << "\n";
156 cout << setw(55) << left << "Sample Standard Deviation" << "= " << Sd << "\n";
157 cout << setw(55) << left << "Expected True Standard Deviation" << "= " << D << "\n\n";
158 //
159 // Now we can calculate and output some stats:
160 //
161 // test-statistic:
162 double t_stat = (N - 1) * (Sd / D) * (Sd / D);
163 cout << setw(55) << left << "Test Statistic" << "= " << t_stat << "\n";
164 //
165 // Finally define our distribution, and get the probability:
166 //
167 chi_squared dist(N - 1);
168 double p = cdf(dist, t_stat);
169 cout << setw(55) << left << "CDF of test statistic: " << "= "
170 << setprecision(3) << scientific << p << "\n";
171 double ucv = quantile(complement(dist, alpha));
172 double ucv2 = quantile(complement(dist, alpha / 2));
173 double lcv = quantile(dist, alpha);
174 double lcv2 = quantile(dist, alpha / 2);
175 cout << setw(55) << left << "Upper Critical Value at alpha: " << "= "
176 << setprecision(3) << scientific << ucv << "\n";
177 cout << setw(55) << left << "Upper Critical Value at alpha/2: " << "= "
178 << setprecision(3) << scientific << ucv2 << "\n";
179 cout << setw(55) << left << "Lower Critical Value at alpha: " << "= "
180 << setprecision(3) << scientific << lcv << "\n";
181 cout << setw(55) << left << "Lower Critical Value at alpha/2: " << "= "
182 << setprecision(3) << scientific << lcv2 << "\n\n";
183 //
184 // Finally print out results of alternative hypothesis:
185 //
186 cout << setw(55) << left <<
187 "Results for Alternative Hypothesis and alpha" << "= "
188 << setprecision(4) << fixed << alpha << "\n\n";
189 cout << "Alternative Hypothesis Conclusion\n";
190 cout << "Standard Deviation != " << setprecision(3) << fixed << D << " ";
191 if((ucv2 < t_stat) || (lcv2 > t_stat))
192 cout << "NOT REJECTED\n";
193 else
194 cout << "REJECTED\n";
195 cout << "Standard Deviation < " << setprecision(3) << fixed << D << " ";
196 if(lcv > t_stat)
197 cout << "NOT REJECTED\n";
198 else
199 cout << "REJECTED\n";
200 cout << "Standard Deviation > " << setprecision(3) << fixed << D << " ";
201 if(ucv < t_stat)
202 cout << "NOT REJECTED\n";
203 else
204 cout << "REJECTED\n";
205 cout << endl << endl;
206 } // void chi_squared_test
207
chi_squared_sample_sized(double diff,double variance)208 void chi_squared_sample_sized(
209 double diff, // difference from variance to detect
210 double variance) // true variance
211 {
212 using namespace std;
213 // using boost::math;
214 using boost::math::chi_squared;
215 using boost::math::quantile;
216 using boost::math::complement;
217 using boost::math::cdf;
218
219 try
220 {
221 cout << // Print out general info:
222 "_____________________________________________________________\n"
223 "Estimated sample sizes required for various confidence levels\n"
224 "_____________________________________________________________\n\n";
225 cout << setprecision(5);
226 cout << setw(40) << left << "True Variance" << "= " << variance << "\n";
227 cout << setw(40) << left << "Difference to detect" << "= " << diff << "\n";
228 //
229 // Define a table of significance levels:
230 //
231 double alpha[] = { 0.5, 0.33333333333333333333333, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
232 //
233 // Print table header:
234 //
235 cout << "\n\n"
236 "_______________________________________________________________\n"
237 "Confidence Estimated Estimated\n"
238 " Value (%) Sample Size Sample Size\n"
239 " (lower one- (upper one-\n"
240 " sided test) sided test)\n"
241 "_______________________________________________________________\n";
242 //
243 // Now print out the data for the table rows.
244 //
245 for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
246 {
247 // Confidence value:
248 cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]);
249 // Calculate df for a lower single-sided test:
250 double df = chi_squared::find_degrees_of_freedom(
251 -diff, alpha[i], alpha[i], variance);
252 // Convert to integral sample size (df is a floating point value in this implementation):
253 double size = ceil(df) + 1;
254 // Print size:
255 cout << fixed << setprecision(0) << setw(16) << right << size;
256 // Calculate df for an upper single-sided test:
257 df = chi_squared::find_degrees_of_freedom(
258 diff, alpha[i], alpha[i], variance);
259 // Convert to integral sample size:
260 size = ceil(df) + 1;
261 // Print size:
262 cout << fixed << setprecision(0) << setw(16) << right << size << endl;
263 }
264 cout << endl;
265 }
266 catch(const std::exception& e)
267 { // Always useful to include try & catch blocks because default policies
268 // are to throw exceptions on arguments that cause errors like underflow, overflow.
269 // Lacking try & catch blocks, the program will abort without a message below,
270 // which may give some helpful clues as to the cause of the exception.
271 std::cout <<
272 "\n""Message from thrown exception was:\n " << e.what() << std::endl;
273 ++error_result;
274 }
275 } // chi_squared_sample_sized
276
main()277 int main()
278 {
279 // Run tests for Gear data
280 // see http://www.itl.nist.gov/div898/handbook/eda/section3/eda3581.htm
281 // Tests measurements of gear diameter.
282 //
283 confidence_limits_on_std_deviation(0.6278908E-02, 100);
284 chi_squared_test(0.6278908E-02, 0.1, 100, 0.05);
285 chi_squared_sample_sized(0.1 - 0.6278908E-02, 0.1);
286 //
287 // Run tests for silicon wafer fabrication data.
288 // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc23.htm
289 // A supplier of 100 ohm.cm silicon wafers claims that his fabrication
290 // process can produce wafers with sufficient consistency so that the
291 // standard deviation of resistivity for the lot does not exceed
292 // 10 ohm.cm. A sample of N = 10 wafers taken from the lot has a
293 // standard deviation of 13.97 ohm.cm
294 //
295 confidence_limits_on_std_deviation(13.97, 10);
296 chi_squared_test(13.97, 10.0, 10, 0.05);
297 chi_squared_sample_sized(13.97 * 13.97 - 100, 100);
298 chi_squared_sample_sized(55, 100);
299 chi_squared_sample_sized(1, 100);
300
301 // List confidence interval multipliers for standard deviation
302 // for a range of numbers of observations from 2 to a million,
303 // and for a few alpha values, 0.1, 0.05, 0.01 for condfidences 90, 95, 99 %
304 confidence_limits_on_std_deviation_alpha(1., 0.1);
305 confidence_limits_on_std_deviation_alpha(1., 0.05);
306 confidence_limits_on_std_deviation_alpha(1., 0.01);
307
308 return error_result;
309 }
310
311 /*
312
313 ________________________________________________
314 2-Sided Confidence Limits For Standard Deviation
315 ________________________________________________
316 Number of Observations = 100
317 Standard Deviation = 0.006278908
318 _____________________________________________
319 Confidence Lower Upper
320 Value (%) Limit Limit
321 _____________________________________________
322 50.000 0.00601 0.00662
323 75.000 0.00582 0.00685
324 90.000 0.00563 0.00712
325 95.000 0.00551 0.00729
326 99.000 0.00530 0.00766
327 99.900 0.00507 0.00812
328 99.990 0.00489 0.00855
329 99.999 0.00474 0.00895
330 ______________________________________________
331 Chi Squared test for sample standard deviation
332 ______________________________________________
333 Number of Observations = 100
334 Sample Standard Deviation = 0.00628
335 Expected True Standard Deviation = 0.10000
336 Test Statistic = 0.39030
337 CDF of test statistic: = 1.438e-099
338 Upper Critical Value at alpha: = 1.232e+002
339 Upper Critical Value at alpha/2: = 1.284e+002
340 Lower Critical Value at alpha: = 7.705e+001
341 Lower Critical Value at alpha/2: = 7.336e+001
342 Results for Alternative Hypothesis and alpha = 0.0500
343 Alternative Hypothesis Conclusion
344 Standard Deviation != 0.100 NOT REJECTED
345 Standard Deviation < 0.100 NOT REJECTED
346 Standard Deviation > 0.100 REJECTED
347 _____________________________________________________________
348 Estimated sample sizes required for various confidence levels
349 _____________________________________________________________
350 True Variance = 0.10000
351 Difference to detect = 0.09372
352 _______________________________________________________________
353 Confidence Estimated Estimated
354 Value (%) Sample Size Sample Size
355 (lower one- (upper one-
356 sided test) sided test)
357 _______________________________________________________________
358 50.000 2 2
359 66.667 2 5
360 75.000 2 10
361 90.000 4 32
362 95.000 5 52
363 99.000 8 102
364 99.900 13 178
365 99.990 18 257
366 99.999 23 337
367 ________________________________________________
368 2-Sided Confidence Limits For Standard Deviation
369 ________________________________________________
370 Number of Observations = 10
371 Standard Deviation = 13.9700000
372 _____________________________________________
373 Confidence Lower Upper
374 Value (%) Limit Limit
375 _____________________________________________
376 50.000 12.41880 17.25579
377 75.000 11.23084 19.74131
378 90.000 10.18898 22.98341
379 95.000 9.60906 25.50377
380 99.000 8.62898 31.81825
381 99.900 7.69466 42.51593
382 99.990 7.04085 55.93352
383 99.999 6.54517 73.00132
384 ______________________________________________
385 Chi Squared test for sample standard deviation
386 ______________________________________________
387 Number of Observations = 10
388 Sample Standard Deviation = 13.97000
389 Expected True Standard Deviation = 10.00000
390 Test Statistic = 17.56448
391 CDF of test statistic: = 9.594e-001
392 Upper Critical Value at alpha: = 1.692e+001
393 Upper Critical Value at alpha/2: = 1.902e+001
394 Lower Critical Value at alpha: = 3.325e+000
395 Lower Critical Value at alpha/2: = 2.700e+000
396 Results for Alternative Hypothesis and alpha = 0.0500
397 Alternative Hypothesis Conclusion
398 Standard Deviation != 10.000 REJECTED
399 Standard Deviation < 10.000 REJECTED
400 Standard Deviation > 10.000 NOT REJECTED
401 _____________________________________________________________
402 Estimated sample sizes required for various confidence levels
403 _____________________________________________________________
404 True Variance = 100.00000
405 Difference to detect = 95.16090
406 _______________________________________________________________
407 Confidence Estimated Estimated
408 Value (%) Sample Size Sample Size
409 (lower one- (upper one-
410 sided test) sided test)
411 _______________________________________________________________
412 50.000 2 2
413 66.667 2 5
414 75.000 2 10
415 90.000 4 32
416 95.000 5 51
417 99.000 7 99
418 99.900 11 174
419 99.990 15 251
420 99.999 20 330
421 _____________________________________________________________
422 Estimated sample sizes required for various confidence levels
423 _____________________________________________________________
424 True Variance = 100.00000
425 Difference to detect = 55.00000
426 _______________________________________________________________
427 Confidence Estimated Estimated
428 Value (%) Sample Size Sample Size
429 (lower one- (upper one-
430 sided test) sided test)
431 _______________________________________________________________
432 50.000 2 2
433 66.667 4 10
434 75.000 8 21
435 90.000 23 71
436 95.000 36 115
437 99.000 71 228
438 99.900 123 401
439 99.990 177 580
440 99.999 232 762
441 _____________________________________________________________
442 Estimated sample sizes required for various confidence levels
443 _____________________________________________________________
444 True Variance = 100.00000
445 Difference to detect = 1.00000
446 _______________________________________________________________
447 Confidence Estimated Estimated
448 Value (%) Sample Size Sample Size
449 (lower one- (upper one-
450 sided test) sided test)
451 _______________________________________________________________
452 50.000 2 2
453 66.667 14696 14993
454 75.000 36033 36761
455 90.000 130079 132707
456 95.000 214283 218612
457 99.000 428628 437287
458 99.900 756333 771612
459 99.990 1095435 1117564
460 99.999 1440608 1469711
461 ________________________________________________
462 2-Sided Confidence Limits For Standard Deviation
463 ________________________________________________
464 Confidence level (two-sided) = 0.1000000
465 Standard Deviation = 1.0000000
466 _____________________________________________
467 Observations Lower Upper
468 Limit Limit
469 _____________________________________________
470 2 0.5102 15.9472
471 3 0.5778 4.4154
472 4 0.6196 2.9200
473 5 0.6493 2.3724
474 6 0.6720 2.0893
475 7 0.6903 1.9154
476 8 0.7054 1.7972
477 9 0.7183 1.7110
478 10 0.7293 1.6452
479 15 0.7688 1.4597
480 20 0.7939 1.3704
481 30 0.8255 1.2797
482 40 0.8454 1.2320
483 50 0.8594 1.2017
484 60 0.8701 1.1805
485 100 0.8963 1.1336
486 120 0.9045 1.1203
487 1000 0.9646 1.0383
488 10000 0.9885 1.0118
489 50000 0.9948 1.0052
490 100000 0.9963 1.0037
491 1000000 0.9988 1.0012
492 ________________________________________________
493 2-Sided Confidence Limits For Standard Deviation
494 ________________________________________________
495 Confidence level (two-sided) = 0.0500000
496 Standard Deviation = 1.0000000
497 _____________________________________________
498 Observations Lower Upper
499 Limit Limit
500 _____________________________________________
501 2 0.4461 31.9102
502 3 0.5207 6.2847
503 4 0.5665 3.7285
504 5 0.5991 2.8736
505 6 0.6242 2.4526
506 7 0.6444 2.2021
507 8 0.6612 2.0353
508 9 0.6755 1.9158
509 10 0.6878 1.8256
510 15 0.7321 1.5771
511 20 0.7605 1.4606
512 30 0.7964 1.3443
513 40 0.8192 1.2840
514 50 0.8353 1.2461
515 60 0.8476 1.2197
516 100 0.8780 1.1617
517 120 0.8875 1.1454
518 1000 0.9580 1.0459
519 10000 0.9863 1.0141
520 50000 0.9938 1.0062
521 100000 0.9956 1.0044
522 1000000 0.9986 1.0014
523 ________________________________________________
524 2-Sided Confidence Limits For Standard Deviation
525 ________________________________________________
526 Confidence level (two-sided) = 0.0100000
527 Standard Deviation = 1.0000000
528 _____________________________________________
529 Observations Lower Upper
530 Limit Limit
531 _____________________________________________
532 2 0.3562 159.5759
533 3 0.4344 14.1244
534 4 0.4834 6.4675
535 5 0.5188 4.3960
536 6 0.5464 3.4848
537 7 0.5688 2.9798
538 8 0.5875 2.6601
539 9 0.6036 2.4394
540 10 0.6177 2.2776
541 15 0.6686 1.8536
542 20 0.7018 1.6662
543 30 0.7444 1.4867
544 40 0.7718 1.3966
545 50 0.7914 1.3410
546 60 0.8065 1.3026
547 100 0.8440 1.2200
548 120 0.8558 1.1973
549 1000 0.9453 1.0609
550 10000 0.9821 1.0185
551 50000 0.9919 1.0082
552 100000 0.9943 1.0058
553 1000000 0.9982 1.0018
554 */
555
556