|
|
Boost-Commit : |
From: pbristow_at_[hidden]
Date: 2007-09-23 13:29:18
Author: pbristow
Date: 2007-09-23 13:29:17 EDT (Sun, 23 Sep 2007)
New Revision: 39487
URL: http://svn.boost.org/trac/boost/changeset/39487
Log:
added code to make a list of CI for an alpha for a range of numbers of observations.
Text files modified:
sandbox/math_toolkit/libs/math/example/chi_square_std_dev_test.cpp | 247 +++++++++++++++++++++++++++++----------
1 files changed, 183 insertions(+), 64 deletions(-)
Modified: sandbox/math_toolkit/libs/math/example/chi_square_std_dev_test.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/example/chi_square_std_dev_test.cpp (original)
+++ sandbox/math_toolkit/libs/math/example/chi_square_std_dev_test.cpp 2007-09-23 13:29:17 EDT (Sun, 23 Sep 2007)
@@ -7,14 +7,18 @@
// 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
@@ -27,8 +31,10 @@
// 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 @@
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 @@
"_____________________________________________\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 @@
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 @@
// 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 @@
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 @@
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<double>:
-//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 @@
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
*/
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk