| // 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> |
| |
| 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 occations. |
| // 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; |
| } |
| } // 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 condfidences 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 0; |
| } |
| |
| /* |
| |
| ________________________________________________ |
| 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 |
| */ |
| |
