chi_squared_examples.qbk

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

[section:cs_eg Chi Squared Distribution Examples]

[section:chi_sq_intervals Confidence Intervals on the Standard Deviation]

Once you have calculated the standard deviation for your data, a legitimate
question to ask is "How reliable is the calculated standard deviation?".
For this situation the Chi Squared distribution can be used to calculate
confidence intervals for the standard deviation.

The full example code & sample output is in
[@../../../example/chi_square_std_dev_test.cpp chi_square_std_deviation_test.cpp].

We'll begin by defining the procedure that will calculate and print out the
confidence intervals:

   void confidence_limits_on_std_deviation(
        double Sd,    // Sample Standard Deviation
        unsigned N)   // Sample size
   {
   
We'll begin by printing out some general information:

   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";

and then define a table of significance levels for which we'll calculate
intervals:

   double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };

The distribution we'll need to calculate the confidence intervals is a
Chi Squared distribution, with N-1 degrees of freedom:

   chi_squared dist(N - 1);

For each value of alpha, the formula for the confidence interval is given by:

[equation chi_squ_tut1]

Where [equation chi_squ_tut2] is the upper critical value, and 
[equation chi_squ_tut3] is the lower critical value of the
Chi Squared distribution.

In code we begin by printing out a table header:

   cout << "\n\n"
           "_____________________________________________\n"
           "Confidence          Lower          Upper\n"
           " Value (%)          Limit          Limit\n"
           "_____________________________________________\n";

and then loop over the values of alpha and calculate the intervals
for each: remember that the lower critical value is the same as the
quantile, and the upper critical value is the same as the quantile
from the complement of the probability:

   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;

To see some example output we'll use the
[@http://www.itl.nist.gov/div898/handbook/eda/section3/eda3581.htm 
gear data] from the __handbook.
The data represents measurements of gear diameter from a manufacturing
process.

[pre'''
________________________________________________
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
''']

So at the 95% confidence level we conclude that the standard deviation
is between 0.00551 and 0.00729.

[h4 Confidence intervals as a function of the number of observations]

Similarly, we can also list the confidence intervals for the standard deviation
for the common confidence levels 95%, for increasing numbers of observations.

The standard deviation used to compute these values is unity,
so the limits listed are *multipliers* for any particular standard deviation.
For example, given a standard deviation of 0.0062789 as in the example 
above; for 100 observations the multiplier is 0.8780
giving the lower confidence limit of 0.8780 * 0.006728 = 0.00551.

[pre'''
____________________________________________________
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
''']

With just 2 observations the limits are from *0.445* up to to *31.9*,
so the standard deviation might be about *half*
the observed value up to *30 times* the observed value!

Estimating a standard deviation with just a handful of values leaves a very great uncertainty,
especially the upper limit.
Note especially how far the upper limit is skewed from the most likely standard deviation.

Even for 10 observations, normally considered a reasonable number,
the range is still from 0.69 to 1.8, about a range of 0.7 to 2,
and is still highly skewed with an upper limit *twice* the median.

When we have 1000 observations, the estimate of the standard deviation is starting to look convincing,
with a range from 0.95 to 1.05 - now near symmetrical, but still about + or - 5%.

Only when we have 10000 or more repeated observations can we start to be reasonably confident
(provided we are sure that other factors like drift are not creeping in).

For 10000 observations, the interval is 0.99 to 1.1 - finally a really convincing + or -1% confidence.

[endsect][/section:chi_sq_intervals Confidence Intervals on the Standard Deviation]

[section:chi_sq_test Chi-Square Test for the Standard Deviation]

We use this test to determine whether the standard deviation of a sample
differs from a specified value.  Typically this occurs in process change
situations where we wish to compare the standard deviation of a new
process to an established one.

The code for this example is contained in 
[@../../../example/chi_square_std_dev_test.cpp chi_square_std_dev_test.cpp], and
we'll begin by defining the procedure that will print out the test
statistics:

   void chi_squared_test(
       double Sd,     // Sample std deviation
       double D,      // True std deviation
       unsigned N,    // Sample size
       double alpha)  // Significance level
   {
   
The procedure begins by printing a summary of the input data:

   using namespace std;
   using namespace boost::math;

   // 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";

The test statistic (T) is simply the ratio of the sample and "true" standard 
deviations squared, multiplied by the number of degrees of freedom (the 
sample size less one):
 
   double t_stat = (N - 1) * (Sd / D) * (Sd / D);
   cout << setw(55) << left << "Test Statistic" << "=  " << t_stat << "\n";

The distribution we need to use, is a Chi Squared distribution with N-1
degrees of freedom:

   chi_squared dist(N - 1);

The various hypothesis that can be tested are summarised in the following table:

[table
[[Hypothesis][Test]]
[[The null-hypothesis: there is no difference in standard deviation from the specified value]
    [ Reject if T < [chi][super 2][sub (1-alpha/2; N-1)] or T > [chi][super 2][sub (alpha/2; N-1)] ]]
[[The alternative hypothesis: there is a difference in standard deviation from the specified value]
    [ Reject if [chi][super 2][sub (1-alpha/2; N-1)] >= T  >= [chi][super 2][sub (alpha/2; N-1)] ]]
[[The alternative hypothesis: the standard deviation is less than the specified value]
    [ Reject if [chi][super 2][sub (1-alpha; N-1)] <= T ]]
[[The alternative hypothesis: the standard deviation is greater than the specified value]
    [ Reject if [chi][super 2][sub (alpha; N-1)] >= T ]]
]

Where [chi][super 2][sub (alpha; N-1)] is the upper critical value of the 
Chi Squared distribution, and [chi][super 2][sub (1-alpha; N-1)] is the
lower critical value.  

Recall that the lower critical value is the same
as the quantile, and the upper critical value is the same as the quantile
from the complement of the probability, that gives us the following code
to calculate the critical values:

   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";