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