// Copyright John Maddock 2006
// Copyright Paul A. Bristow 2007, 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)

#ifdef _MSC_VER
#  pragma warning(disable: 4512) // assignment operator could not be generated.
#  pragma warning(disable: 4510) // default constructor could not be generated.
#  pragma warning(disable: 4610) // can never be instantiated - user defined constructor required.
#endif

#include <boost/math/distributions/students_t.hpp>

// avoid "using namespace std;" and "using namespace boost::math;"
// to avoid potential ambiguity with names in std random.
#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;

void confidence_limits_on_mean(double Sm, double Sd, unsigned Sn)
{
   //
   // Sm = Sample Mean.
   // Sd = Sample Standard Deviation.
   // Sn = Sample Size.
   //
   // Calculate confidence intervals for the mean.
   // 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 mean
   // 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 mean.
   // The interval computed from a given sample either
   // contains the true mean or it does not.
   // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm

   using boost::math::students_t;

   // Print out general info:
   cout <<
      "__________________________________\n"
      "2-Sided Confidence Limits For Mean\n"
      "__________________________________\n\n";
   cout << setprecision(7);
   cout << setw(40) << left << "Number of Observations" << "=  " << Sn << "\n";
   cout << setw(40) << left << "Mean" << "=  " << Sm << "\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:
   //
   students_t dist(Sn - 1);
   //
   // Print table header:
   //
   cout << "\n\n"
           "_______________________________________________________________\n"
           "Confidence       T           Interval          Lower          Upper\n"
           " Value (%)     Value          Width            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 T:
      double T = quantile(complement(dist, alpha[i] / 2));
      // Print T:
      cout << fixed << setprecision(3) << setw(10) << right << T;
      // Calculate width of interval (one sided):
      double w = T * Sd / sqrt(double(Sn));
      // Print width:
      if(w < 0.01)
         cout << scientific << setprecision(3) << setw(17) << right << w;
      else
         cout << fixed << setprecision(3) << setw(17) << right << w;
      // Print Limits:
      cout << fixed << setprecision(5) << setw(15) << right << Sm - w;
      cout << fixed << setprecision(5) << setw(15) << right << Sm + w << endl;
   }
   cout << endl;
} // void confidence_limits_on_mean

void single_sample_t_test(double M, double Sm, double Sd, unsigned Sn, double alpha)
{
   //
   // M = true mean.
   // Sm = Sample Mean.
   // Sd = Sample Standard Deviation.
   // Sn = Sample Size.
   // alpha = Significance Level.
   //
   // A Students t test applied to a single set of data.
   // We are testing the null hypothesis that the true
   // mean of the sample is M, 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/eda352.htm
   
   using boost::math::students_t;

   // Print header:
   cout <<
      "__________________________________\n"
      "Student t test for a single sample\n"
      "__________________________________\n\n";
   cout << setprecision(5);
   cout << setw(55) << left << "Number of Observations" << "=  " << Sn << "\n";
   cout << setw(55) << left << "Sample Mean" << "=  " << Sm << "\n";
   cout << setw(55) << left << "Sample Standard Deviation" << "=  " << Sd << "\n";
   cout << setw(55) << left << "Expected True Mean" << "=  " << M << "\n\n";
   //
   // Now we can calculate and output some stats:
   //
   // Difference in means:
   double diff = Sm - M;
   cout << setw(55) << left << "Sample Mean - Expected Test Mean" << "=  " << diff << "\n";
   // Degrees of freedom:
   unsigned v = Sn - 1;
   cout << setw(55) << left << "Degrees of Freedom" << "=  " << v << "\n";
   // t-statistic:
   double t_stat = diff * sqrt(double(Sn)) / Sd;
   cout << setw(55) << left << "T Statistic" << "=  " << t_stat << "\n";
   //
   // Finally define our distribution, and get the probability:
   //
   students_t dist(v);
   double q = cdf(complement(dist, fabs(t_stat)));
   cout << setw(55) << left << "Probability that difference is due to chance" << "=  "
      << setprecision(3) << scientific << 2 * q << "\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 << "Mean != " << setprecision(3) << fixed << M << "            ";
   if(q < alpha / 2)
      cout << "NOT REJECTED\n";
   else
      cout << "REJECTED\n";
   cout << "Mean  < " << setprecision(3) << fixed << M << "            ";
   if(cdf(dist, t_stat) < alpha)
      cout << "NOT REJECTED\n";
   else
      cout << "REJECTED\n";
   cout << "Mean  > " << setprecision(3) << fixed << M << "            ";
   if(cdf(complement(dist, t_stat)) < alpha)
      cout << "NOT REJECTED\n";
   else
      cout << "REJECTED\n";
   cout << endl << endl;
} // void single_sample_t_test(

void single_sample_find_df(double M, double Sm, double Sd)
{
   //
   // M = true mean.
   // Sm = Sample Mean.
   // Sd = Sample Standard Deviation.
   //
 
   using boost::math::students_t;

   // Print out general info:
   cout <<
      "_____________________________________________________________\n"
      "Estimated sample sizes required for various confidence levels\n"
      "_____________________________________________________________\n\n";
   cout << setprecision(5);
   cout << setw(40) << left << "True Mean" << "=  " << M << "\n";
   cout << setw(40) << left << "Sample Mean" << "=  " << Sm << "\n";
   cout << setw(40) << left << "Sample Standard Deviation" << "=  " << Sd << "\n";
   //
   // Define a table of significance intervals:
   //
   double alpha[] = { 0.5, 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"
           "              (one sided test)    (two sided test)\n"
           "_______________________________________________________________\n";
   //
   // Now print out the data for the table rows.
   //
   for(unsigned i = 1; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
   {
      // Confidence value:
      cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]);
      // calculate df for single sided test:
      double df = students_t::find_degrees_of_freedom(
         fabs(M - Sm), alpha[i], alpha[i], Sd);
      // convert to sample size, always one more than the degrees of freedom:
      double size = ceil(df) + 1;
      // Print size:
      cout << fixed << setprecision(0) << setw(16) << right << size;
      // calculate df for two sided test:
      df = students_t::find_degrees_of_freedom(
         fabs(M - Sm), alpha[i]/2, alpha[i], Sd);
      // convert to sample size:
      size = ceil(df) + 1;
      // Print size:
      cout << fixed << setprecision(0) << setw(16) << right << size << endl;
   }
   cout << endl;
} // void single_sample_find_df

int main()
{
   //
   // Run tests for Heat Flow Meter data
   // see http://www.itl.nist.gov/div898/handbook/eda/section4/eda428.htm
   // The data was collected while calibrating a heat flow meter
   // against a known value.
   //
   confidence_limits_on_mean(9.261460, 0.2278881e-01, 195);
   single_sample_t_test(5, 9.261460, 0.2278881e-01, 195, 0.05);
   single_sample_find_df(5, 9.261460, 0.2278881e-01);

   //
   // Data for this example from:
   // P.K.Hou, O. W. Lau & M.C. Wong, Analyst (1983) vol. 108, p 64.
   // from Statistics for Analytical Chemistry, 3rd ed. (1994), pp 54-55
   // J. C. Miller and J. N. Miller, Ellis Horwood ISBN 0 13 0309907
   //
   // Determination of mercury by cold-vapour atomic absorption,
   // the following values were obtained fusing a trusted
   // Standard Reference Material containing 38.9% mercury,
   // which we assume is correct or 'true'.
   //
   confidence_limits_on_mean(37.8, 0.964365, 3);
   // 95% test:
   single_sample_t_test(38.9, 37.8, 0.964365, 3, 0.05);
   // 90% test:
   single_sample_t_test(38.9, 37.8, 0.964365, 3, 0.1);
   // parameter estimate:
   single_sample_find_df(38.9, 37.8, 0.964365);

   return 0;
} // int main()

/*

Output:

------ Rebuild All started: Project: students_t_single_sample, Configuration: Release Win32 ------
  students_t_single_sample.cpp
  Generating code
  Finished generating code
  students_t_single_sample.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\students_t_single_sample.exe
  __________________________________
  2-Sided Confidence Limits For Mean
  __________________________________
  
  Number of Observations                  =  195
  Mean                                    =  9.26146
  Standard Deviation                      =  0.02278881
  
  
  _______________________________________________________________
  Confidence       T           Interval          Lower          Upper
   Value (%)     Value          Width            Limit          Limit
  _______________________________________________________________
      50.000     0.676       1.103e-003        9.26036        9.26256
      75.000     1.154       1.883e-003        9.25958        9.26334
      90.000     1.653       2.697e-003        9.25876        9.26416
      95.000     1.972       3.219e-003        9.25824        9.26468
      99.000     2.601       4.245e-003        9.25721        9.26571
      99.900     3.341       5.453e-003        9.25601        9.26691
      99.990     3.973       6.484e-003        9.25498        9.26794
      99.999     4.537       7.404e-003        9.25406        9.26886
  
  __________________________________
  Student t test for a single sample
  __________________________________
  
  Number of Observations                                 =  195
  Sample Mean                                            =  9.26146
  Sample Standard Deviation                              =  0.02279
  Expected True Mean                                     =  5.00000
  
  Sample Mean - Expected Test Mean                       =  4.26146
  Degrees of Freedom                                     =  194
  T Statistic                                            =  2611.28380
  Probability that difference is due to chance           =  0.000e+000
  
  Results for Alternative Hypothesis and alpha           =  0.0500
  
  Alternative Hypothesis     Conclusion
  Mean != 5.000            NOT REJECTED
  Mean  < 5.000            REJECTED
  Mean  > 5.000            NOT REJECTED
  
  
  _____________________________________________________________
  Estimated sample sizes required for various confidence levels
  _____________________________________________________________
  
  True Mean                               =  5.00000
  Sample Mean                             =  9.26146
  Sample Standard Deviation               =  0.02279
  
  
  _______________________________________________________________
  Confidence       Estimated          Estimated
   Value (%)      Sample Size        Sample Size
                (one sided test)    (two sided test)
  _______________________________________________________________
      75.000               2               2
      90.000               2               2
      95.000               2               2
      99.000               2               2
      99.900               3               3
      99.990               3               3
      99.999               4               4
  
  __________________________________
  2-Sided Confidence Limits For Mean
  __________________________________
  
  Number of Observations                  =  3
  Mean                                    =  37.8000000
  Standard Deviation                      =  0.9643650
  
  
  _______________________________________________________________
  Confidence       T           Interval          Lower          Upper
   Value (%)     Value          Width            Limit          Limit
  _______________________________________________________________
      50.000     0.816            0.455       37.34539       38.25461
      75.000     1.604            0.893       36.90717       38.69283
      90.000     2.920            1.626       36.17422       39.42578
      95.000     4.303            2.396       35.40438       40.19562
      99.000     9.925            5.526       32.27408       43.32592
      99.900    31.599           17.594       20.20639       55.39361
      99.990    99.992           55.673      -17.87346       93.47346
      99.999   316.225          176.067     -138.26683      213.86683
  
  __________________________________
  Student t test for a single sample
  __________________________________
  
  Number of Observations                                 =  3
  Sample Mean                                            =  37.80000
  Sample Standard Deviation                              =  0.96437
  Expected True Mean                                     =  38.90000
  
  Sample Mean - Expected Test Mean                       =  -1.10000
  Degrees of Freedom                                     =  2
  T Statistic                                            =  -1.97566
  Probability that difference is due to chance           =  1.869e-001
  
  Results for Alternative Hypothesis and alpha           =  0.0500
  
  Alternative Hypothesis     Conclusion
  Mean != 38.900            REJECTED
  Mean  < 38.900            REJECTED
  Mean  > 38.900            REJECTED
  
  
  __________________________________
  Student t test for a single sample
  __________________________________
  
  Number of Observations                                 =  3
  Sample Mean                                            =  37.80000
  Sample Standard Deviation                              =  0.96437
  Expected True Mean                                     =  38.90000
  
  Sample Mean - Expected Test Mean                       =  -1.10000
  Degrees of Freedom                                     =  2
  T Statistic                                            =  -1.97566
  Probability that difference is due to chance           =  1.869e-001
  
  Results for Alternative Hypothesis and alpha           =  0.1000
  
  Alternative Hypothesis     Conclusion
  Mean != 38.900            REJECTED
  Mean  < 38.900            NOT REJECTED
  Mean  > 38.900            REJECTED
  
  
  _____________________________________________________________
  Estimated sample sizes required for various confidence levels
  _____________________________________________________________
  
  True Mean                               =  38.90000
  Sample Mean                             =  37.80000
  Sample Standard Deviation               =  0.96437
  
  
  _______________________________________________________________
  Confidence       Estimated          Estimated
   Value (%)      Sample Size        Sample Size
                (one sided test)    (two sided test)
  _______________________________________________________________
      75.000               3               4
      90.000               7               9
      95.000              11              13
      99.000              20              22
      99.900              35              37
      99.990              50              53
      99.999              66              68
  

*/