libs/math/example/chi_square_std_dev_test.cpp
// 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 */