From d262d6d4c34df2b13bfdcb28dc902573ddd4908d Mon Sep 17 00:00:00 2001 From: "Paul A. Bristow" Date: Sun, 23 Sep 2007 17:29:17 +0000 Subject: [PATCH] added code to make a list of CI for an alpha for a range of numbers of observations. [SVN r39487] --- example/chi_square_std_dev_test.cpp | 247 +++++++++++++++++++++------- 1 file changed, 183 insertions(+), 64 deletions(-) diff --git a/example/chi_square_std_dev_test.cpp b/example/chi_square_std_dev_test.cpp index a06063fda..005c861a0 100644 --- a/example/chi_square_std_dev_test.cpp +++ b/example/chi_square_std_dev_test.cpp @@ -7,14 +7,18 @@ // or copy at http://www.boost.org/LICENSE_1_0.txt) #include +using std::cout; using std::endl; +using std::left; using std::fixed; using std::right; using std::scientific; #include +using std::setw; +using std::setprecision; #include + 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 @@ -27,8 +31,10 @@ void confidence_limits_on_std_deviation( // contains the true standard deviation or it does not. // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda358.htm - using namespace std; - using namespace boost::math; + // using namespace boost::math; + using boost::math::chi_squared; + using boost::math::quantile; + using boost::math::complement; // Print out general info: cout << @@ -40,11 +46,9 @@ void confidence_limits_on_std_deviation( 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: @@ -56,7 +60,6 @@ void confidence_limits_on_std_deviation( "_____________________________________________\n"; // // Now print out the data for the table rows. - // for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i) { // Confidence value: @@ -69,7 +72,57 @@ void confidence_limits_on_std_deviation( 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 @@ -85,8 +138,11 @@ void chi_squared_test( // that any difference is not down to chance. // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda358.htm // - using namespace std; - using namespace boost::math; + // using namespace boost::math; + using boost::math::chi_squared; + using boost::math::quantile; + using boost::math::complement; + using boost::math::cdf; // Print header: cout << @@ -145,21 +201,23 @@ void chi_squared_test( 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 namespace boost::math; + // using boost::math; + using boost::math::chi_squared; + using boost::math::quantile; + using boost::math::complement; + using boost::math::cdf; try { - - // Print out general info: - cout << - "_____________________________________________________________\n" + cout << // Print out general info: + "_____________________________________________________________\n" "Estimated sample sizes required for various confidence levels\n" "_____________________________________________________________\n\n"; cout << setprecision(5); @@ -211,16 +269,10 @@ void chi_squared_sample_sized( std::cout << "\n""Message from thrown exception was:\n " << e.what() << std::endl; } -} - -//Message from thrown exception was: for alpha 0.5 -// Error in function boost::math::tools::bracket_and_solve_root: -//Unable to bracket root, last nearest value was 1.2446030555722283e-058 - +} // 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. @@ -242,21 +294,24 @@ int main() 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; } /* -Output: - ________________________________________________ 2-Sided Confidence Limits For Standard Deviation ________________________________________________ - Number of Observations = 100 Standard Deviation = 0.006278908 - - _____________________________________________ Confidence Lower Upper Value (%) Limit Limit @@ -269,45 +324,36 @@ _____________________________________________ 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 + (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 @@ -315,15 +361,11 @@ _______________________________________________________________ 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 @@ -336,45 +378,36 @@ _____________________________________________ 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 + (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 @@ -382,22 +415,19 @@ _______________________________________________________________ 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 + (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 @@ -405,22 +435,19 @@ _______________________________________________________________ 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 + (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 @@ -428,6 +455,98 @@ _______________________________________________________________ 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 */