6e8145a076
Change the heuristics used when the initial guess is really bad. Change example to return a fail-code when an error occurs.
556 lines
24 KiB
C++
556 lines
24 KiB
C++
// 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 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;
|
|
++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 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 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
|
|
*/
|
|
|