Changeset 80723


Ignore:
Timestamp:
Sep 27, 2012, 4:19:08 PM (6 years ago)
Author:
Paul A. Bristow
Message:

Added trac #7402 resolution

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/libs/math/doc/sf_and_dist/distributions/students_t_examples.qbk

    r51455 r80723  
    66Let's say you have a sample mean, you may wish to know what confidence intervals
    77you can place on that mean.  Colloquially: "I want an interval that I can be
    8 P% sure contains the true mean".  (On a technical point, note that 
    9 the interval either contains the true mean or it does not: the 
    10 meaning of the confidence level is subtly 
    11 different from this colloquialism.  More background information can be found on the 
     8P% sure contains the true mean".  (On a technical point, note that
     9the interval either contains the true mean or it does not: the
     10meaning of the confidence level is subtly
     11different from this colloquialism.  More background information can be found on the
    1212[@http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm NIST site]).
    1313
     
    1616[equation dist_tutorial4]
    1717
    18 Where, ['Y[sub s]] is the sample mean, /s/ is the sample standard deviation, 
    19 /N/ is the sample size, /[alpha]/ is the desired significance level and 
     18Where, ['Y[sub s]] is the sample mean, /s/ is the sample standard deviation,
     19/N/ is the sample size, /[alpha]/ is the desired significance level and
    2020['t[sub ([alpha]/2,N-1)]] is the upper critical value of the Students-t
    2121distribution with /N-1/ degrees of freedom.
     
    2828The confidence level of the test is defined as 1 - [alpha], and often expressed
    2929as a percentage.  So for example a significance level of 0.05, is equivalent
    30 to a 95% confidence level.  Refer to 
    31 [@http://www.itl.nist.gov/div898/handbook/prc/section1/prc14.htm 
     30to a 95% confidence level.  Refer to
     31[@http://www.itl.nist.gov/div898/handbook/prc/section1/prc14.htm
    3232"What are confidence intervals?"] in __handbook for more information.
    3333] [/ Note]
     
    5050[@../../../example/students_t_single_sample.cpp students_t_single_sample.cpp].
    5151
    52 We'll begin by defining a procedure to calculate intervals for 
     52We'll begin by defining a procedure to calculate intervals for
    5353various confidence levels; the procedure will print these out
    5454as a table:
     
    6060   // Bring everything into global namespace for ease of use:
    6161   using namespace boost::math;
    62    using namespace std;   
     62   using namespace std;
    6363
    6464   void confidence_limits_on_mean(
     
    7171
    7272      // Print out general info:
    73       cout << 
     73      cout <<
    7474         "__________________________________\n"
    7575         "2-Sided Confidence Limits For Mean\n"
     
    8383
    8484      double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
    85      
     85
    8686Note that these are the complements of the confidence/probability levels: 0.5, 0.75, 0.9 .. 0.99999).
    8787
     
    9898
    9999   double T = quantile(complement(dist, alpha[i] / 2));
    100    
     100
    101101Note that alpha was divided by two, since we'll be calculating
    102102both the upper and lower bounds: had we been interested in a single
    103103sided interval then we would have omitted this step.
    104    
     104
    105105Now to complete the picture, we'll get the (one-sided) width of the
    106106interval from the t-statistic
     
    116116Let's take a look at some sample output, first using the
    117117[@http://www.itl.nist.gov/div898/handbook/eda/section4/eda428.htm
    118 Heat flow data] from the NIST site.  The data set was collected 
    119 by Bob Zarr of NIST in January, 1990 from a heat flow meter 
     118Heat flow data] from the NIST site.  The data set was collected
     119by Bob Zarr of NIST in January, 1990 from a heat flow meter
    120120calibration and stability analysis.
    121121The corresponding dataplot
    122 output for this test can be found in 
    123 [@http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm 
     122output for this test can be found in
     123[@http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm
    124124section 3.5.2] of the __handbook.
    125125
     
    147147       99.999     4.537       7.404e-003        9.25406        9.26886
    148148''']
    149    
     149
    150150As you can see the large sample size (195) and small standard deviation (0.023)
    151151have combined to give very small intervals, indeed we can be
     
    156156and from Statistics for Analytical Chemistry, 3rd ed. (1994), pp 54-55
    157157J. C. Miller and J. N. Miller, Ellis Horwood ISBN 0 13 0309907.]
    158 The values result from the determination of mercury by cold-vapour 
     158The values result from the determination of mercury by cold-vapour
    159159atomic absorption.
    160160
     
    191191[section:tut_mean_test Testing a sample mean for difference from a "true" mean]
    192192
    193 When calibrating or comparing a scientific instrument or measurement method of some kind, 
     193When calibrating or comparing a scientific instrument or measurement method of some kind,
    194194we want to be answer the question "Does an observed sample mean differ from the
    195195"true" mean in any significant way?".  If it does, then we have evidence of
    196196a systematic difference.  This question can be answered with a Students-t test:
    197 more information can be found 
    198 [@http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm 
     197more information can be found
     198[@http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm
    199199on the NIST site].
    200200
     
    214214but statisticans eschew this to avoid implying that there is positive evidence of 'no difference'.
    215215'Not-rejected' here means there is *no evidence* of difference, but there still might well be a difference.
    216 For example, see [@http://en.wikipedia.org/wiki/Argument_from_ignorance argument from ignorance] and 
     216For example, see [@http://en.wikipedia.org/wiki/Argument_from_ignorance argument from ignorance] and
    217217[@http://www.bmj.com/cgi/content/full/311/7003/485 Absence of evidence does not constitute evidence of absence.]
    218218] [/ note]
     
    225225   // Bring everything into global namespace for ease of use:
    226226   using namespace boost::math;
    227    using namespace std;   
     227   using namespace std;
    228228
    229229   void single_sample_t_test(double M, double Sm, double Sd, unsigned Sn, double alpha)
     
    235235      // Sn = Sample Size.
    236236      // alpha = Significance Level.
    237      
     237
    238238Most of the procedure is pretty-printing, so let's just focus on the
    239239calculation, we begin by calculating the t-statistic:
     
    246246   double t_stat = diff * sqrt(double(Sn)) / Sd;
    247247
    248 Finally calculate the probability from the t-statistic. If we're interested 
     248Finally calculate the probability from the t-statistic. If we're interested
    249249in simply whether there is a difference (either less or greater) or not,
    250250we don't care about the sign of the t-statistic,
     
    254254   students_t dist(v);
    255255   double q = cdf(complement(dist, fabs(t_stat)));
    256    
     256
    257257The procedure then prints out the results of the various tests
    258 that can be done, these 
     258that can be done, these
    259259can be summarised in the following table:
    260260
     
    262262[[Hypothesis][Test]]
    263263[[The Null-hypothesis: there is
    264 *no difference* in means] 
     264*no difference* in means]
    265265[Reject if complement of CDF for |t| < significance level / 2:
    266266
     
    268268
    269269[[The Alternative-hypothesis: there
    270 *is difference* in means] 
     270*is difference* in means]
    271271[Reject if complement of CDF for |t| > significance level / 2:
    272272
     
    274274
    275275[[The Alternative-hypothesis: the sample mean *is less* than
    276 the true mean.] 
     276the true mean.]
    277277[Reject if CDF of t > significance level:
    278278
     
    280280
    281281[[The Alternative-hypothesis: the sample mean *is greater* than
    282 the true mean.] 
     282the true mean.]
    283283[Reject if complement of CDF of t > significance level:
    284284
     
    290290and against `alpha` for a one-sided test]
    291291
    292 Now that we have all the parts in place, let's take a look at some 
     292Now that we have all the parts in place, let's take a look at some
    293293sample output, first using the
    294294[@http://www.itl.nist.gov/div898/handbook/eda/section4/eda428.htm
    295 Heat flow data] from the NIST site.  The data set was collected 
    296 by Bob Zarr of NIST in January, 1990 from a heat flow meter 
     295Heat flow data] from the NIST site.  The data set was collected
     296by Bob Zarr of NIST in January, 1990 from a heat flow meter
    297297calibration and stability analysis.  The corresponding dataplot
    298 output for this test can be found in 
    299 [@http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm 
     298output for this test can be found in
     299[@http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm
    300300section 3.5.2] of the __handbook.
    301301
     
    335335and from Statistics for Analytical Chemistry, 3rd ed. (1994), pp 54-55
    336336J. C. Miller and J. N. Miller, Ellis Horwood ISBN 0 13 0309907.]
    337 The values result from the determination of mercury by cold-vapour 
     337The values result from the determination of mercury by cold-vapour
    338338atomic absorption.
    339339
     
    364364in the location of the true mean.  So even though there appears to be a difference
    365365between the sample mean and the expected true mean, we conclude that there
    366 is no significant difference, and are unable to reject the null hypothesis. 
    367 However, if we were to lower the bar for acceptance down to alpha = 0.1 
     366is no significant difference, and are unable to reject the null hypothesis.
     367However, if we were to lower the bar for acceptance down to alpha = 0.1
    368368(a 90% confidence level) we see a different output:
    369369
     
    407407can provide this information.
    408408
    409 This section is based on the example code in 
     409This section is based on the example code in
    410410[@../../../example/students_t_single_sample.cpp students_t_single_sample.cpp]
    411411and we begin by defining a procedure that will print out a table of
     
    418418   // Bring everything into global namespace for ease of use:
    419419   using namespace boost::math;
    420    using namespace std;   
     420   using namespace std;
    421421
    422422   void single_sample_find_df(
     
    448448The first argument is the difference between the means that you
    449449wish to be able to detect, here it's the absolute value of the
    450 difference between the sample mean, and the true mean. 
     450difference between the sample mean, and the true mean.
    451451
    452452Then come two probability values: alpha and beta.  Alpha is the
     
    454454in fact true.  Beta is the maximum acceptable risk of failing to reject
    455455the null-hypothesis when in fact it is false.
    456 Also note that for a two-sided test, alpha must be divided by 2. 
     456Also note that for a two-sided test, alpha must be divided by 2.
    457457
    458458The final parameter of the function is the standard deviation of the sample.
     
    488488and from Statistics for Analytical Chemistry, 3rd ed. (1994), pp 54-55
    489489J. C. Miller and J. N. Miller, Ellis Horwood ISBN 0 13 0309907.]
    490 The values result from the determination of mercury by cold-vapour 
    491 atomic absorption. 
     490The values result from the determination of mercury by cold-vapour
     491atomic absorption.
    492492
    493493Only three measurements were made, and the Students-t test above
     
    519519''']
    520520
    521 So in this case, many more measurements would have had to be made, 
     521So in this case, many more measurements would have had to be made,
    522522for example at the 95% level, 14 measurements in total for a two-sided test.
    523523
     
    529529determining whether a new process or treatment is better than an old one.
    530530
    531 In this example, we'll be using the 
     531In this example, we'll be using the
    532532[@http://www.itl.nist.gov/div898/handbook/eda/section3/eda3531.htm
    533 Car Mileage sample data] from the 
     533Car Mileage sample data] from the
    534534[@http://www.itl.nist.gov NIST website].  The data compares
    535535miles per gallon of US cars with miles per gallon of Japanese cars.
    536536
    537 The sample code is in 
     537The sample code is in
    538538[@../../../example/students_t_two_samples.cpp students_t_two_samples.cpp].
    539539
     
    542542If the standard deviations are assumed to be equal, then the calculation
    543543of the t-statistic is greatly simplified, so we'll examine that case first.
    544 In real life we should verify whether this assumption is valid with a 
     544In real life we should verify whether this assumption is valid with a
    545545Chi-Squared test for equal variances.
    546546
     
    565565           double alpha)     // alpha = Significance Level.
    566566   {
    567      
     567
    568568
    569569Our procedure will begin by calculating the t-statistic, assuming
     
    572572[equation dist_tutorial1]
    573573
    574 where Sp is the "pooled" standard deviation of the two samples, 
     574where Sp is the "pooled" standard deviation of the two samples,
    575575and /v/ is the number of degrees of freedom of the two combined
    576576samples.  We can now write the code to calculate the t-statistic:
     
    581581   // Pooled variance:
    582582   double sp = sqrt(((Sn1-1) * Sd1 * Sd1 + (Sn2-1) * Sd2 * Sd2) / v);
    583    cout << setw(55) << left << "Pooled Standard Deviation" << "=  " << v << "\n";
     583   cout << setw(55) << left << "Pooled Standard Deviation" << "=  " << sp << "\n";
    584584   // t-statistic:
    585585   double t_stat = (Sm1 - Sm2) / (sp * sqrt(1.0 / Sn1 + 1.0 / Sn2));
    586586   cout << setw(55) << left << "T Statistic" << "=  " << t_stat << "\n";
    587    
     587
    588588The next step is to define our distribution object, and calculate the
    589589complement of the probability:
     
    591591   students_t dist(v);
    592592   double q = cdf(complement(dist, fabs(t_stat)));
    593    cout << setw(55) << left << "Probability that difference is due to chance" << "=  " 
     593   cout << setw(55) << left << "Probability that difference is due to chance" << "=  "
    594594      << setprecision(3) << scientific << 2 * q << "\n\n";
    595595
     
    603603[[Hypothesis][Test]]
    604604[[The Null-hypothesis: there is
    605 *no difference* in means] 
     605*no difference* in means]
    606606[Reject if complement of CDF for |t| < significance level / 2:
    607607
     
    609609
    610610[[The Alternative-hypothesis: there is a
    611 *difference* in means] 
     611*difference* in means]
    612612[Reject if complement of CDF for |t| > significance level / 2:
    613613
     
    615615
    616616[[The Alternative-hypothesis: Sample 1 Mean is *less* than
    617 Sample 2 Mean.] 
     617Sample 2 Mean.]
    618618[Reject if CDF of t > significance level:
    619619
     
    621621
    622622[[The Alternative-hypothesis: Sample 1 Mean is *greater* than
    623 Sample 2 Mean.] 
     623Sample 2 Mean.]
    624624
    625625[Reject if complement of CDF of t > significance level:
     
    634634skip over that, and take a look at the sample output for alpha=0.05
    635635(a 95% probability level).  For comparison the dataplot output
    636 for the same data is in 
     636for the same data is in
    637637[@http://www.itl.nist.gov/div898/handbook/eda/section3/eda353.htm
    638638section 1.3.5.3] of the __handbook.
     
    644644
    645645   Number of Observations (Sample 1)                      =  249
    646    Sample 1 Mean                                          =  20.14458
    647    Sample 1 Standard Deviation                            =  6.41470
     646   Sample 1 Mean                                          =  20.145
     647   Sample 1 Standard Deviation                            =  6.4147
    648648   Number of Observations (Sample 2)                      =  79
    649    Sample 2 Mean                                          =  30.48101
    650    Sample 2 Standard Deviation                            =  6.10771
    651    Degrees of Freedom                                     =  326.00000
    652    Pooled Standard Deviation                              =  326.00000
    653    T Statistic                                            =  -12.62059
     649   Sample 2 Mean                                          =  30.481
     650   Sample 2 Standard Deviation                            =  6.1077
     651   Degrees of Freedom                                     =  326
     652   Pooled Standard Deviation                              =  6.3426
     653   T Statistic                                            =  -12.621
    654654   Probability that difference is due to chance           =  5.273e-030
    655655
     
    667667The tests on the alternative hypothesis show that we must
    668668also reject the hypothesis that Sample 1 Mean is
    669 greater than that for Sample 2: in this case Sample 1 represents the 
    670 miles per gallon for Japanese cars, and Sample 2 the miles per gallon for 
     669greater than that for Sample 2: in this case Sample 1 represents the
     670miles per gallon for Japanese cars, and Sample 2 the miles per gallon for
    671671US cars, so we conclude that Japanese cars are on average more
    672672fuel efficient.
     
    678678[equation dist_tutorial2]
    679679
    680 And for the combined degrees of freedom we use the 
     680And for the combined degrees of freedom we use the
    681681[@http://en.wikipedia.org/wiki/Welch-Satterthwaite_equation Welch-Satterthwaite]
    682682approximation:
     
    692692an integer value: this may be necessary if you are relying on lookup tables,
    693693but since our code fully supports non-integer degrees of freedom there is no
    694 need to truncate in this case.  Also note that when the degrees of freedom 
     694need to truncate in this case.  Also note that when the degrees of freedom
    695695is small then the Welch-Satterthwaite approximation may be a significant
    696696source of error.]
     
    754754Instead we calculate the difference between before and after measurements
    755755in each patient, and calculate the mean and standard deviation of the differences.
    756 To test whether a significant change has taken place, we can then test 
    757 the null-hypothesis that the true mean is zero using the same procedure 
     756To test whether a significant change has taken place, we can then test
     757the null-hypothesis that the true mean is zero using the same procedure
    758758we used in the single sample cases previously discussed.
    759759
    760760That means we can:
    761761
    762 * [link math_toolkit.dist.stat_tut.weg.st_eg.tut_mean_intervals Calculate confidence intervals of the mean]. 
    763 If the endpoints of the interval differ in sign then we are unable to reject 
     762* [link math_toolkit.dist.stat_tut.weg.st_eg.tut_mean_intervals Calculate confidence intervals of the mean].
     763If the endpoints of the interval differ in sign then we are unable to reject
    764764the null-hypothesis that there is no change.
    765765* [link math_toolkit.dist.stat_tut.weg.st_eg.tut_mean_test Test whether the true mean is zero]. If the
     
    773773[endsect][/section:st_eg Student's t]
    774774
    775 [/ 
    776   Copyright 2006 John Maddock and Paul A. Bristow.
     775[/
     776  Copyright 2006, 2012 John Maddock and Paul A. Bristow.
    777777  Distributed under the Boost Software License, Version 1.0.
    778778  (See accompanying file LICENSE_1_0.txt or copy at
Note: See TracChangeset for help on using the changeset viewer.