123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500 |
- [section:cs_eg Chi Squared Distribution Examples]
- [section:chi_sq_intervals Confidence Intervals on the Standard Deviation]
- Once you have calculated the standard deviation for your data, a legitimate
- question to ask is "How reliable is the calculated standard deviation?".
- For this situation the Chi Squared distribution can be used to calculate
- confidence intervals for the standard deviation.
- The full example code & sample output is in
- [@../../example/chi_square_std_dev_test.cpp chi_square_std_dev_test.cpp].
- We'll begin by defining the procedure that will calculate and print out the
- confidence intervals:
- void confidence_limits_on_std_deviation(
- double Sd, // Sample Standard Deviation
- unsigned N) // Sample size
- {
- We'll begin by printing out some general information:
- 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";
- and then define a table of significance levels for which we'll calculate
- intervals:
- double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
- The distribution we'll need to calculate the confidence intervals is a
- Chi Squared distribution, with N-1 degrees of freedom:
- chi_squared dist(N - 1);
- For each value of alpha, the formula for the confidence interval is given by:
- [equation chi_squ_tut1]
- Where [equation chi_squ_tut2] is the upper critical value, and
- [equation chi_squ_tut3] is the lower critical value of the
- Chi Squared distribution.
- In code we begin by printing out a table header:
- cout << "\n\n"
- "_____________________________________________\n"
- "Confidence Lower Upper\n"
- " Value (%) Limit Limit\n"
- "_____________________________________________\n";
- and then loop over the values of alpha and calculate the intervals
- for each: remember that the lower critical value is the same as the
- quantile, and the upper critical value is the same as the quantile
- from the complement of the probability:
- 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;
- To see some example output we'll use the
- [@http://www.itl.nist.gov/div898/handbook/eda/section3/eda3581.htm
- gear data] from the __handbook.
- The data represents measurements of gear diameter from a manufacturing
- process.
- [pre'''
- ________________________________________________
- 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
- ''']
- So at the 95% confidence level we conclude that the standard deviation
- is between 0.00551 and 0.00729.
- [h4 Confidence intervals as a function of the number of observations]
- Similarly, we can also list the confidence intervals for the standard deviation
- for the common confidence levels 95%, for increasing numbers of observations.
- The standard deviation used to compute these values is unity,
- so the limits listed are *multipliers* for any particular standard deviation.
- For example, given a standard deviation of 0.0062789 as in the example
- above; for 100 observations the multiplier is 0.8780
- giving the lower confidence limit of 0.8780 * 0.006728 = 0.00551.
- [pre'''
- ____________________________________________________
- 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
- ''']
- With just 2 observations the limits are from *0.445* up to to *31.9*,
- so the standard deviation might be about *half*
- the observed value up to [*30 times] the observed value!
- Estimating a standard deviation with just a handful of values leaves a very great uncertainty,
- especially the upper limit.
- Note especially how far the upper limit is skewed from the most likely standard deviation.
- Even for 10 observations, normally considered a reasonable number,
- the range is still from 0.69 to 1.8, about a range of 0.7 to 2,
- and is still highly skewed with an upper limit *twice* the median.
- When we have 1000 observations, the estimate of the standard deviation is starting to look convincing,
- with a range from 0.95 to 1.05 - now near symmetrical, but still about + or - 5%.
- Only when we have 10000 or more repeated observations can we start to be reasonably confident
- (provided we are sure that other factors like drift are not creeping in).
- For 10000 observations, the interval is 0.99 to 1.1 - finally a really convincing + or -1% confidence.
- [endsect] [/section:chi_sq_intervals Confidence Intervals on the Standard Deviation]
- [section:chi_sq_test Chi-Square Test for the Standard Deviation]
- We use this test to determine whether the standard deviation of a sample
- differs from a specified value. Typically this occurs in process change
- situations where we wish to compare the standard deviation of a new
- process to an established one.
- The code for this example is contained in
- [@../../example/chi_square_std_dev_test.cpp chi_square_std_dev_test.cpp], and
- we'll begin by defining the procedure that will print out the test
- statistics:
- void chi_squared_test(
- double Sd, // Sample std deviation
- double D, // True std deviation
- unsigned N, // Sample size
- double alpha) // Significance level
- {
- The procedure begins by printing a summary of the input data:
- using namespace std;
- using namespace boost::math;
- // 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";
- The test statistic (T) is simply the ratio of the sample and "true" standard
- deviations squared, multiplied by the number of degrees of freedom (the
- sample size less one):
- double t_stat = (N - 1) * (Sd / D) * (Sd / D);
- cout << setw(55) << left << "Test Statistic" << "= " << t_stat << "\n";
- The distribution we need to use, is a Chi Squared distribution with N-1
- degrees of freedom:
- chi_squared dist(N - 1);
- The various hypothesis that can be tested are summarised in the following table:
- [table
- [[Hypothesis][Test]]
- [[The null-hypothesis: there is no difference in standard deviation from the specified value]
- [ Reject if T < [chi][super 2][sub (1-alpha/2; N-1)] or T > [chi][super 2][sub (alpha/2; N-1)] ]]
- [[The alternative hypothesis: there is a difference in standard deviation from the specified value]
- [ Reject if [chi][super 2][sub (1-alpha/2; N-1)] >= T >= [chi][super 2][sub (alpha/2; N-1)] ]]
- [[The alternative hypothesis: the standard deviation is less than the specified value]
- [ Reject if [chi][super 2][sub (1-alpha; N-1)] <= T ]]
- [[The alternative hypothesis: the standard deviation is greater than the specified value]
- [ Reject if [chi][super 2][sub (alpha; N-1)] >= T ]]
- ]
- Where [chi][super 2][sub (alpha; N-1)] is the upper critical value of the
- Chi Squared distribution, and [chi][super 2][sub (1-alpha; N-1)] is the
- lower critical value.
- Recall that the lower critical value is the same
- as the quantile, and the upper critical value is the same as the quantile
- from the complement of the probability, that gives us the following code
- to calculate the critical values:
- 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";
- Now that we have the critical values, we can compare these to our test
- statistic, and print out the result of each hypothesis and test:
- 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 << "ACCEPTED\n";
- else
- cout << "REJECTED\n";
- cout << "Standard Deviation < " << setprecision(3) << fixed << D << " ";
- if(lcv > t_stat)
- cout << "ACCEPTED\n";
- else
- cout << "REJECTED\n";
- cout << "Standard Deviation > " << setprecision(3) << fixed << D << " ";
- if(ucv < t_stat)
- cout << "ACCEPTED\n";
- else
- cout << "REJECTED\n";
- cout << endl << endl;
- To see some example output we'll use the
- [@http://www.itl.nist.gov/div898/handbook/eda/section3/eda3581.htm
- gear data] from the __handbook.
- The data represents measurements of gear diameter from a manufacturing
- process. The program output is deliberately designed to mirror
- the DATAPLOT output shown in the
- [@http://www.itl.nist.gov/div898/handbook/eda/section3/eda358.htm
- NIST Handbook Example].
- [pre'''
- ______________________________________________
- 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 ACCEPTED
- Standard Deviation < 0.100 ACCEPTED
- Standard Deviation > 0.100 REJECTED
- ]
- In this case we are testing whether the sample standard deviation is 0.1,
- and the null-hypothesis is rejected, so we conclude that the standard
- deviation ['is not] 0.1.
- For an alternative example, consider the
- [@http://www.itl.nist.gov/div898/handbook/prc/section2/prc23.htm
- silicon wafer data] again from the __handbook.
- In this scenario 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, and the question
- we ask ourselves is "Is the suppliers claim correct?".
- The program output now looks like this:
- [pre'''
- ______________________________________________
- 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 ACCEPTED
- ]
- In this case, our null-hypothesis is that the standard deviation of
- the sample is less than 10: this hypothesis is rejected in the analysis
- above, and so we reject the manufacturers claim.
- [endsect] [/section:chi_sq_test Chi-Square Test for the Standard Deviation]
- [section:chi_sq_size Estimating the Required Sample Sizes for a Chi-Square Test for the Standard Deviation]
- Suppose we conduct a Chi Squared test for standard deviation and the result
- is borderline, a legitimate question to ask is "How large would the sample size
- have to be in order to produce a definitive result?"
- The class template [link math_toolkit.dist_ref.dists.chi_squared_dist
- chi_squared_distribution] has a static method
- `find_degrees_of_freedom` that will calculate this value for
- some acceptable risk of type I failure /alpha/, type II failure
- /beta/, and difference from the standard deviation /diff/. Please
- note that the method used works on variance, and not standard deviation
- as is usual for the Chi Squared Test.
- The code for this example is located in
- [@../../example/chi_square_std_dev_test.cpp chi_square_std_dev_test.cpp].
- We begin by defining a procedure to print out the sample sizes required
- for various risk levels:
- void chi_squared_sample_sized(
- double diff, // difference from variance to detect
- double variance) // true variance
- {
- The procedure begins by printing out the input data:
- using namespace std;
- using namespace boost::math;
- // Print out general info:
- cout <<
- "_____________________________________________________________\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";
- And defines a table of significance levels for which we'll calculate sample sizes:
- double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
- For each value of alpha we can calculate two sample sizes: one where the
- sample variance is less than the true value by /diff/ and one
- where it is greater than the true value by /diff/. Thanks to the
- asymmetric nature of the Chi Squared distribution these two values will
- not be the same, the difference in their calculation differs only in the
- sign of /diff/ that's passed to `find_degrees_of_freedom`. Finally
- in this example we'll simply things, and let risk level /beta/ be the
- same as /alpha/:
- 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 sample size:
- 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 sample size:
- size = ceil(df) + 1;
- // Print size:
- cout << fixed << setprecision(0) << setw(16) << right << size << endl;
- }
- cout << endl;
- For some example output, consider the
- [@http://www.itl.nist.gov/div898/handbook/prc/section2/prc23.htm
- silicon wafer data] from the __handbook.
- In this scenario 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, and the question
- we ask ourselves is "How large would our sample have to be to reliably
- detect this difference?".
- To use our procedure above, we have to convert the
- standard deviations to variance (square them),
- after which the program output looks like this:
- [pre'''
- _____________________________________________________________
- 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
- 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'''
- ]
- In this case we are interested in a upper single sided test.
- So for example, if the maximum acceptable risk of falsely rejecting
- the null-hypothesis is 0.05 (Type I error), and the maximum acceptable
- risk of failing to reject the null-hypothesis is also 0.05
- (Type II error), we estimate that we would need a sample size of 51.
- [endsect] [/section:chi_sq_size Estimating the Required Sample Sizes for a Chi-Square Test for the Standard Deviation]
- [endsect] [/section:cs_eg Chi Squared Distribution]
- [/
- Copyright 2006, 2013 John Maddock and Paul A. Bristow.
- Distributed under 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).
- ]
|