// Copyright Paul A. Bristow 2006.
// Copyright John Maddock 2006.

// 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)

// test_students_t.cpp

// http://en.wikipedia.org/wiki/Student%27s_t_distribution
// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3664.htm

// Basic sanity test for Student's t probability (quantile) (0. < p < 1).
// and Student's t probability Quantile (0. < p < 1).

#include <boost/test/test_exec_monitor.hpp> // Boost.Test
#include <boost/test/floating_point_comparison.hpp>

#include <boost/math/concepts/real_concept.hpp> // for real_concept
#include <boost/math/distributions/students_t.hpp>
    using boost::math::students_t_distribution;
#include <boost/math/tools/test.hpp> // for real_concept

#include <iostream>
   using std::cout;
   using std::endl;
   using std::setprecision;
#include <limits>
  using std::numeric_limits;

template <class RealType>
RealType naive_pdf(RealType v, RealType t)
{
   // Calculate the pdf of the students t in a deliberately
   // naive way, using equation (5) from
   // http://mathworld.wolfram.com/Studentst-Distribution.html
   // This is equivalent to, but a different method
   // to the one in the actual implementation, so can be used as
   // a very basic sanity check.  However some published values
   // would be nice....

   using namespace std;  // for ADL
   using boost::math::beta;

   //return pow(v / (v + t*t), (1+v) / 2) / (sqrt(v) * beta(v/2, RealType(0.5f)));
   RealType result = boost::math::tgamma_ratio((v+1)/2, v/2);
   result /= sqrt(v * boost::math::constants::pi<RealType>());
   result /= pow(1 + t*t/v, (v+1)/2);
   return result;
}

template <class RealType>
void test_spots(RealType)
{
  // Basic sanity checks

   RealType tolerance = static_cast<RealType>(1e-4); // 1e-6 (as %)
   // Some tests only pass at 1e-5 because probability value is less accurate,
   // a digit in 6th decimal place, although calculated using
   // a t-distribution generator (claimed 6 decimal digits) at
  // http://faculty.vassar.edu/lowry/VassarStats.html
   // http://faculty.vassar.edu/lowry/tsamp.html
   // df = 5, +/-t = 2.0, 1-tailed = 0.050970, 2-tailed = 0.101939

   cout << "Tolerance for type " << typeid(RealType).name()  << " is " << tolerance << " %" << endl;

   // http://en.wikipedia.org/wiki/Student%27s_t_distribution#Table_of_selected_values
  // Using tabulated value of t = 3.182 for 0.975, 3 df, one-sided.

   // http://www.mth.kcl.ac.uk/~shaww/web_page/papers/Tdistribution06.pdf refers to:

   // A lookup table of quantiles of the RealType distribution
  // for 1 to 25 in steps of 0.1 is provided in CSV form at:
  // www.mth.kcl.ac.uk/~shaww/web_page/papers/Tsupp/tquantiles.csv
   // gives accurate t of -3.1824463052837 and 3 degrees of freedom.
   // Values below are from this source, saved as tquantiles.xls.
   // DF are across the columns, probabilities down the rows
   // and the t- values (quantiles) are shown.
   // These values are probably accurate to nearly 64-bit double
  // (perhaps 14 decimal digits).

   BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(2),       // degrees_of_freedom
         static_cast<RealType>(-6.96455673428326)),  // t
         static_cast<RealType>(0.01),                // probability.
         tolerance); // %

   BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(5),       // degrees_of_freedom
         static_cast<RealType>(-3.36492999890721)),  // t
         static_cast<RealType>(0.01),                // probability.
         tolerance);

   BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(1),      // degrees_of_freedom
         static_cast<RealType>(-31830.988607907)),  // t
         static_cast<RealType>(0.00001),            // probability.
         tolerance);

   BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(25.),    // degrees_of_freedom
         static_cast<RealType>(-5.2410429995425)),  // t
         static_cast<RealType>(0.00001),            // probability.
         tolerance);

   BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(1),   // degrees_of_freedom
         static_cast<RealType>(-63661.97723)),   // t
         static_cast<RealType>(0.000005),        // probability.
         tolerance);

    BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(5.),  // degrees_of_freedom
         static_cast<RealType>(-17.89686614)),   // t
         static_cast<RealType>(0.000005),        // probability.
         tolerance);

    BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(25.),  // degrees_of_freedom
         static_cast<RealType>(-5.510848412)),    // t
         static_cast<RealType>(0.000005),         // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(10.),  // degrees_of_freedom
         static_cast<RealType>(-1.812461123)),    // t
         static_cast<RealType>(0.05),             // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(10),  // degrees_of_freedom
         static_cast<RealType>(1.812461123)),    // t
         static_cast<RealType>(0.95),            // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         complement(
            students_t_distribution<RealType>(10),  // degrees_of_freedom
            static_cast<RealType>(1.812461123))),    // t
         static_cast<RealType>(0.05),            // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(10),  // degrees_of_freedom
         static_cast<RealType>(9.751995491)),    // t
         static_cast<RealType>(0.999999),        // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(10.),  // degrees_of_freedom - for ALL degrees_of_freedom!
         static_cast<RealType>(0.)),              // t
         static_cast<RealType>(0.5),              // probability.
         tolerance);


   // Student's t Inverse function tests.
  // Special cases

  BOOST_CHECK_THROW(boost::math::quantile(
         students_t_distribution<RealType>(1.),  // degrees_of_freedom (ignored).
         static_cast<RealType>(0)), std::overflow_error); // t == -infinity.

  BOOST_CHECK_THROW(boost::math::quantile(
         students_t_distribution<RealType>(1.),  // degrees_of_freedom (ignored).
         static_cast<RealType>(1)), std::overflow_error); // t == +infinity.

  BOOST_CHECK_EQUAL(boost::math::quantile(
         students_t_distribution<RealType>(1.),  // degrees_of_freedom (ignored).
         static_cast<RealType>(0.5)),  //  probability == half - special case.
         static_cast<RealType>(0)); // t == zero.

  BOOST_CHECK_EQUAL(boost::math::quantile(
         complement(
            students_t_distribution<RealType>(1.),  // degrees_of_freedom (ignored).
            static_cast<RealType>(0.5))),  //  probability == half - special case.
         static_cast<RealType>(0)); // t == zero.

  BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(1.),  // degrees_of_freedom (ignored).
         static_cast<RealType>(0.5)),  //  probability == half - special case.
         static_cast<RealType>(0), // t == zero.
         tolerance);

   BOOST_CHECK_CLOSE( // Tests of p middling.
      ::boost::math::cdf(
         students_t_distribution<RealType>(5.),  // degrees_of_freedom
         static_cast<RealType>(-0.559429644)),  // t
         static_cast<RealType>(0.3), // probability.
         tolerance);
   BOOST_CHECK_CLOSE(
      ::boost::math::quantile(
         students_t_distribution<RealType>(5.),  // degrees_of_freedom
         static_cast<RealType>(0.3)),  // probability.
         static_cast<RealType>(-0.559429644), // t
         tolerance);

   BOOST_CHECK_CLOSE(
      ::boost::math::quantile(
         complement(
            students_t_distribution<RealType>(5.),  // degrees_of_freedom
            static_cast<RealType>(0.7))),  // probability.
         static_cast<RealType>(-0.559429644), // t
         tolerance);

   BOOST_CHECK_CLOSE( // Tests of p high.
      ::boost::math::cdf(
         students_t_distribution<RealType>(5.),  // degrees_of_freedom
         static_cast<RealType>(1.475884049)),  // t
         static_cast<RealType>(0.9), // probability.
         tolerance);
   BOOST_CHECK_CLOSE(
      ::boost::math::quantile(
         students_t_distribution<RealType>(5.),  // degrees_of_freedom
         static_cast<RealType>(0.9)),  // probability.
         static_cast<RealType>(1.475884049), // t
         tolerance);

   BOOST_CHECK_CLOSE( // Tests of p low.
      ::boost::math::cdf(
         students_t_distribution<RealType>(5.),  // degrees_of_freedom
         static_cast<RealType>(-1.475884049)),  // t
         static_cast<RealType>(0.1), // probability.
         tolerance);
   BOOST_CHECK_CLOSE(
      ::boost::math::quantile(
         students_t_distribution<RealType>(5.),  // degrees_of_freedom
         static_cast<RealType>(0.1)),  // probability.
         static_cast<RealType>(-1.475884049), // t
         tolerance);

   BOOST_CHECK_CLOSE(
      ::boost::math::cdf(
         students_t_distribution<RealType>(2.),  // degrees_of_freedom
         static_cast<RealType>(-6.96455673428326)),  // t
         static_cast<RealType>(0.01), // probability.
         tolerance);

   BOOST_CHECK_CLOSE(
      ::boost::math::quantile(
         students_t_distribution<RealType>(2.),  // degrees_of_freedom
         static_cast<RealType>(0.01)),  // probability.
         static_cast<RealType>(-6.96455673428326), // t
         tolerance);

      //
      // Some special tests to exercise the double-precision approximations
      // to the quantile:
      //
      // tolerance is 50 eps expressed as a persent:
      //
      tolerance = boost::math::tools::epsilon<RealType>() * 5000;
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(2.00390625L),                     // degrees_of_freedom.
         static_cast<RealType>(0.5625L)),                                    //  probability.
         static_cast<RealType>(0.178133131573788108465134803511798566L),     // t.
         tolerance);      
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(1L),                              // degrees_of_freedom.
         static_cast<RealType>(0.03125L)),                                   //  probability.
         static_cast<RealType>(-10.1531703876088604621071476634194722L),     // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(1L),                            // degrees_of_freedom.
         static_cast<RealType>(0.875L)),                                   //  probability.
         static_cast<RealType>(2.41421356237309504880168872421390942L),    // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(2L),                              // degrees_of_freedom.
         static_cast<RealType>(0.03125L)),                                   //  probability.
         static_cast<RealType>(-3.81000381000571500952501666878143315L),     // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(2L),                            // degrees_of_freedom.
         static_cast<RealType>(0.875L)),                                   //  probability.
         static_cast<RealType>(1.60356745147454630810732088527854144L),    // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(4L),                              // degrees_of_freedom.
         static_cast<RealType>(0.03125L)),                                   //  probability.
         static_cast<RealType>(-2.56208431914409044861223047927635034L),     // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(4L),                            // degrees_of_freedom.
         static_cast<RealType>(0.875L)),                                   //  probability.
         static_cast<RealType>(1.34439755550909142430681981315923574L),    // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(6L),                              // degrees_of_freedom.
         static_cast<RealType>(0.03125L)),                                   //  probability.
         static_cast<RealType>(-2.28348667906973065861212495010082952L),     // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(6L),                            // degrees_of_freedom.
         static_cast<RealType>(0.875L)),                                   //  probability.
         static_cast<RealType>(1.27334930914664286821103236660071906L),    // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(8L),                              // degrees_of_freedom.
         static_cast<RealType>(0.03125L)),                                   //  probability.
         static_cast<RealType>(-2.16296475406014719458642055768894376L),     // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(8L),                            // degrees_of_freedom.
         static_cast<RealType>(0.875L)),                                   //  probability.
         static_cast<RealType>(1.24031826078267310637634677726479038L),    // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(10L),                             // degrees_of_freedom.
         static_cast<RealType>(0.03125L)),                                   //  probability.
         static_cast<RealType>(-2.09596136475109350926340169211429572L),     // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(10L),                         // degrees_of_freedom.
         static_cast<RealType>(0.875L)),                                 //  probability.
         static_cast<RealType>(1.2212553950039221407185188573696834L),   // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(2.125L),                          // degrees_of_freedom.
         static_cast<RealType>(0.03125L)),                                   //  probability.
         static_cast<RealType>(-3.62246031671091980110493455859296532L),     // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(2.125L),                        // degrees_of_freedom.
         static_cast<RealType>(0.875L)),                                   //  probability.
         static_cast<RealType>(1.56905270993307293450392958697861969L),    // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(3L),                              // degrees_of_freedom.
         static_cast<RealType>(0.03125L)),                                   //  probability.
         static_cast<RealType>(-2.90004411882995814036141778367917946L),     // t.
         tolerance);
      BOOST_CHECK_CLOSE(boost::math::quantile(
         students_t_distribution<RealType>(3L),                            // degrees_of_freedom.
         static_cast<RealType>(0.875L)),                                   //  probability.
         static_cast<RealType>(1.42262528146180931868169289781115099L),    // t.
         tolerance);

  // Student's t pdf tests.
  // for PDF checks, use 100 eps tolerance expressed as a percent:
   tolerance = boost::math::tools::epsilon<RealType>() * 10000;

   for(unsigned i = 1; i < 20; i += 3)
   {
      for(RealType r = -10; r < 10; r += 0.125)
      {
         //std::cout << "df=" << i << " t=" << r << std::endl;
         BOOST_CHECK_CLOSE(
            boost::math::pdf(
               students_t_distribution<RealType>(static_cast<RealType>(i)),
               r),
            naive_pdf<RealType>(static_cast<RealType>(i), r),
            tolerance);
      }
   }

    RealType tol2 = boost::math::tools::epsilon<RealType>() * 5;
    students_t_distribution<RealType> dist(8);
    RealType x = static_cast<RealType>(0.125);
    using namespace std; // ADL of std names.
    // mean:
    BOOST_CHECK_CLOSE(
       mean(dist)
       , static_cast<RealType>(0), tol2);
    // variance:
    BOOST_CHECK_CLOSE(
       variance(dist)
       , static_cast<RealType>(8.0L / 6.0L), tol2);
    // std deviation:
    BOOST_CHECK_CLOSE(
       standard_deviation(dist)
       , static_cast<RealType>(sqrt(8.0L / 6.0L)), tol2);
    // hazard:
    BOOST_CHECK_CLOSE(
       hazard(dist, x)
       , pdf(dist, x) / cdf(complement(dist, x)), tol2);
    // cumulative hazard:
    BOOST_CHECK_CLOSE(
       chf(dist, x)
       , -log(cdf(complement(dist, x))), tol2);
    // coefficient_of_variation:
    BOOST_CHECK_THROW(
       coefficient_of_variation(dist),
       std::overflow_error);
    // mode:
    BOOST_CHECK_CLOSE(
       mean(dist)
       , static_cast<RealType>(0), tol2);
    // median:
    BOOST_CHECK_CLOSE(
       median(dist)
       , static_cast<RealType>(0), tol2);
    // skewness:
    BOOST_CHECK_CLOSE(
       skewness(dist)
       , static_cast<RealType>(0), tol2);
    // kurtosis:
    BOOST_CHECK_CLOSE(
       kurtosis(dist)
       , static_cast<RealType>(4.5), tol2);
    // kurtosis excess:
    BOOST_CHECK_CLOSE(
       kurtosis_excess(dist)
       , static_cast<RealType>(1.5), tol2);

    // Parameter estimation. These results are close to but
    // not identical to those reported on the NIST website at
    // http://www.itl.nist.gov/div898/handbook/prc/section2/prc222.htm
    // the NIST results appear to be calculated using a normal
    // approximation, which slightly under-estimates the degrees of
    // freedom required, particularly when the result is small.
    //
    BOOST_CHECK_EQUAL(
       ceil(students_t_distribution<RealType>::find_degrees_of_freedom(
         static_cast<RealType>(0.5),
         static_cast<RealType>(0.005),
         static_cast<RealType>(0.01),
         static_cast<RealType>(1.0))),
         99);
    BOOST_CHECK_EQUAL(
       ceil(students_t_distribution<RealType>::find_degrees_of_freedom(
         static_cast<RealType>(1.5),
         static_cast<RealType>(0.005),
         static_cast<RealType>(0.01),
         static_cast<RealType>(1.0))),
         14);
    BOOST_CHECK_EQUAL(
       ceil(students_t_distribution<RealType>::find_degrees_of_freedom(
         static_cast<RealType>(0.5),
         static_cast<RealType>(0.025),
         static_cast<RealType>(0.01),
         static_cast<RealType>(1.0))),
         76);
    BOOST_CHECK_EQUAL(
       ceil(students_t_distribution<RealType>::find_degrees_of_freedom(
         static_cast<RealType>(1.5),
         static_cast<RealType>(0.025),
         static_cast<RealType>(0.01),
         static_cast<RealType>(1.0))),
         11);
    BOOST_CHECK_EQUAL(
       ceil(students_t_distribution<RealType>::find_degrees_of_freedom(
         static_cast<RealType>(0.5),
         static_cast<RealType>(0.05),
         static_cast<RealType>(0.01),
         static_cast<RealType>(1.0))),
         65);
    BOOST_CHECK_EQUAL(
       ceil(students_t_distribution<RealType>::find_degrees_of_freedom(
         static_cast<RealType>(1.5),
         static_cast<RealType>(0.05),
         static_cast<RealType>(0.01),
         static_cast<RealType>(1.0))),
         9);

} // template <class RealType>void test_spots(RealType)

int test_main(int, char* [])
{
  // Check that can construct students_t distribution using the two convenience methods:
  using namespace boost::math;
  students_t myst1(2); // Using typedef
   students_t_distribution<> myst2(2); // Using default RealType double.
   //students_t_distribution<double> myst3(2); // Using explicit RealType double.

    // Basic sanity-check spot values.
   // (Parameter value, arbitrarily zero, only communicates the floating point type).
  test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
  test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  test_spots(0.0L); // Test long double.
#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
  test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
#endif
#else
   std::cout << "<note>The long double tests have been disabled on this platform "
      "either because the long double overloads of the usual math functions are "
      "not available at all, or because they are too inaccurate for these tests "
      "to pass.</note>" << std::cout;
#endif

   return 0;
} // int test_main(int, char* [])

/*

Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_students_t.exe"
Running 1 test case...
Tolerance for type float is 0.0001 %
Tolerance for type double is 0.0001 %
Tolerance for type long double is 0.0001 %
Tolerance for type class boost::math::concepts::real_concept is 0.0001 %
*** No errors detected

*/