1 // students_t_example1.cpp
2
3 // Copyright Paul A. Bristow 2006, 2007.
4
5 // Use, modification and distribution are subject to the
6 // Boost Software License, Version 1.0.
7 // (See accompanying file LICENSE_1_0.txt
8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
9
10 // Example 1 of using Student's t
11
12 // http://en.wikipedia.org/wiki/Student's_t-test says:
13 // The t statistic was invented by William Sealy Gosset
14 // for cheaply monitoring the quality of beer brews.
15 // "Student" was his pen name.
16 // WS Gosset was statistician for Guinness brewery in Dublin, Ireland,
17 // hired due to Claude Guinness's innovative policy of recruiting the
18 // best graduates from Oxford and Cambridge for applying biochemistry
19 // and statistics to Guinness's industrial processes.
20 // Gosset published the t test in Biometrika in 1908,
21 // but was forced to use a pen name by his employer who regarded the fact
22 // that they were using statistics as a trade secret.
23 // In fact, Gosset's identity was unknown not only to fellow statisticians
24 // but to his employer - the company insisted on the pseudonym
25 // so that it could turn a blind eye to the breach of its rules.
26
27 // Data for this example from:
28 // P.K.Hou, O. W. Lau & M.C. Wong, Analyst (1983) vol. 108, p 64.
29 // from Statistics for Analytical Chemistry, 3rd ed. (1994), pp 54-55
30 // J. C. Miller and J. N. Miller, Ellis Horwood ISBN 0 13 0309907
31
32 // Determination of mercury by cold-vapour atomic absorption,
33 // the following values were obtained fusing a trusted
34 // Standard Reference Material containing 38.9% mercury,
35 // which we assume is correct or 'true'.
36 double standard = 38.9;
37
38 const int values = 3;
39 double value[values] = {38.9, 37.4, 37.1};
40
41 // Is there any evidence for systematic error?
42
43 // The Students't distribution function is described at
44 // http://en.wikipedia.org/wiki/Student%27s_t_distribution
45 #include <boost/math/distributions/students_t.hpp>
46 using boost::math::students_t; // Probability of students_t(df, t).
47
48 #include <iostream>
49 using std::cout; using std::endl;
50 #include <iomanip>
51 using std::setprecision;
52 #include <cmath>
53 using std::sqrt;
54
main()55 int main()
56 {
57 cout << "Example 1 using Student's t function. " << endl;
58
59 // Example/test using tabulated value
60 // (deliberately coded as naively as possible).
61
62 // Null hypothesis is that there is no difference (greater or less)
63 // between measured and standard.
64
65 double degrees_of_freedom = values-1; // 3-1 = 2
66 cout << "Measurement 1 = " << value[0] << ", measurement 2 = " << value[1] << ", measurement 3 = " << value[2] << endl;
67 double mean = (value[0] + value[1] + value[2]) / static_cast<double>(values);
68 cout << "Standard = " << standard << ", mean = " << mean << ", (mean - standard) = " << mean - standard << endl;
69 double sd = sqrt(((value[0] - mean) * (value[0] - mean) + (value[1] - mean) * (value[1] - mean) + (value[2] - mean) * (value[2] - mean))/ static_cast<double>(values-1));
70 cout << "Standard deviation = " << sd << endl;
71 if (sd == 0.)
72 {
73 cout << "Measured mean is identical to SRM value," << endl;
74 cout << "so probability of no difference between measured and standard (the 'null hypothesis') is unity." << endl;
75 return 0;
76 }
77
78 double t = (mean - standard) * std::sqrt(static_cast<double>(values)) / sd;
79 cout << "Student's t = " << t << endl;
80 cout.precision(2); // Useful accuracy is only a few decimal digits.
81 cout << "Probability of Student's t is " << cdf(students_t(degrees_of_freedom), std::abs(t)) << endl;
82 // 0.91, is 1 tailed.
83 // So there is insufficient evidence of a difference to meet a 95% (1 in 20) criterion.
84
85 return 0;
86 } // int main()
87
88 /*
89
90 Output is:
91
92 Example 1 using Student's t function.
93 Measurement 1 = 38.9, measurement 2 = 37.4, measurement 3 = 37.1
94 Standard = 38.9, mean = 37.8, (mean - standard) = -1.1
95 Standard deviation = 0.964365
96 Student's t = -1.97566
97 Probability of Student's t is 0.91
98
99 */
100
101
102