test_nc_chi_squared.cpp

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

// test_nc_chi_squared.cpp

// Copyright John Maddock 2008.

// 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:4127 4512)
#endif

#if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
#  define TEST_FLOAT
#  define TEST_DOUBLE
#  define TEST_LDOUBLE
#  define TEST_REAL_CONCEPT
#endif

#include <boost/math/concepts/real_concept.hpp> // for real_concept
#include <boost/math/distributions/non_central_chi_squared.hpp> // for chi_squared_distribution
#include <boost/math/special_functions/cbrt.hpp> // for chi_squared_distribution
#include <boost/test/included/test_exec_monitor.hpp> // for test_main
#include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE

#include "functor.hpp"
#include "handle_test_result.hpp"
#include "test_nccs_hooks.hpp"

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

#define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
   {\
      unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
      BOOST_CHECK_CLOSE(a, b, prec); \
      if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
      {\
         std::cerr << "Failure was at row " << i << std::endl;\
         std::cerr << std::setprecision(35); \
         std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
         std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
      }\
   }

#define BOOST_CHECK_EX(a, i) \
   {\
      unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
      BOOST_CHECK(a); \
      if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
      {\
         std::cerr << "Failure was at row " << i << std::endl;\
         std::cerr << std::setprecision(35); \
         std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
         std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
      }\
   }

void expected_results()
{
   //
   // Define the max and mean errors expected for
   // various compilers and platforms.
   //
   const char* largest_type;
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
   if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
   {
      largest_type = "(long\\s+)?double|real_concept";
   }
   else
   {
      largest_type = "long double|real_concept";
   }
#else
   largest_type = "(long\\s+)?double|real_concept";
#endif

   add_expected_result(
      "[^|]*",                          // compiler
      "[^|]*",                          // stdlib
      "Mac OS",                          // platform
      largest_type,                     // test type(s)
      "[^|]*medium[^|]*",                   // test data group
      "[^|]*", 550, 100);                  // test function
   //
   // Catch all cases come last:
   //
   add_expected_result(
      "[^|]*",                          // compiler
      "[^|]*",                          // stdlib
      "[^|]*",                          // platform
      largest_type,                     // test type(s)
      "[^|]*medium[^|]*",                   // test data group
      "[^|]*", 350, 100);                  // test function
   add_expected_result(
      "[^|]*",                          // compiler
      "[^|]*",                          // stdlib
      "[^|]*",                          // platform
      largest_type,                     // test type(s)
      "[^|]*large[^|]*",                   // test data group
      "[^|]*", 17000, 3000);                  // test function

   //
   // Allow some long double error to creep into
   // the double results:
   //
   add_expected_result(
      "[^|]*",                          // compiler
      "[^|]*",                          // stdlib
      "[^|]*",                          // platform
      "double",                         // test type(s)
      "[^|]*",                   // test data group
      "[^|]*", 3, 2);                  // test function

   //
   // Finish off by printing out the compiler/stdlib/platform names,
   // we do this to make it easier to mark up expected error rates.
   //
   std::cout << "Tests run with " << BOOST_COMPILER << ", " 
      << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
}

template <class RealType>
RealType naive_pdf(RealType v, RealType lam, RealType x)
{
   // Formula direct from 
   // http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
   // with no simplification:
   RealType sum, term, prefix(1);
   RealType eps = boost::math::tools::epsilon<RealType>();
   term = sum = pdf(boost::math::chi_squared_distribution<RealType>(v), x);
   for(int i = 1;; ++i)
   {
      prefix *= lam / (2 * i);
      term = prefix * pdf(boost::math::chi_squared_distribution<RealType>(v + 2 * i), x);
      sum += term;
      if(term / sum < eps)
         break;
   }
   return sum * exp(-lam/2);
}

template <class RealType>
void test_spot(
     RealType df,    // Degrees of freedom
     RealType ncp,   // non-centrality param
     RealType cs,    // Chi Square statistic
     RealType P,     // CDF
     RealType Q,     // Complement of CDF
     RealType tol)   // Test tolerance
{
   boost::math::non_central_chi_squared_distribution<RealType> dist(df, ncp);
   BOOST_CHECK_CLOSE(
      cdf(dist, cs), P, tol);
   try{
      BOOST_CHECK_CLOSE(
         pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), ncp, cs), tol * 50);
   }
   catch(const std::overflow_error&)
   {}
   if((P < 0.99) && (Q < 0.99))
   {
      //
      // We can only check this if P is not too close to 1,
      // so that we can guarentee Q is reasonably free of error:
      //
      BOOST_CHECK_CLOSE(
         cdf(complement(dist, cs)), Q, tol);
      BOOST_CHECK_CLOSE(
            quantile(dist, P), cs, tol * 10);
      BOOST_CHECK_CLOSE(
            quantile(complement(dist, Q)), cs, tol * 10);
      BOOST_CHECK_CLOSE(
         dist.find_degrees_of_freedom(ncp, cs, P), df, tol * 10);
      BOOST_CHECK_CLOSE(
         dist.find_degrees_of_freedom(boost::math::complement(ncp, cs, Q)), df, tol * 10);
      BOOST_CHECK_CLOSE(
         dist.find_non_centrality(df, cs, P), ncp, tol * 10);
      BOOST_CHECK_CLOSE(
         dist.find_non_centrality(boost::math::complement(df, cs, Q)), ncp, tol * 10);
   }
}

template <class RealType> // Any floating-point type RealType.
void test_spots(RealType)
{
   RealType tolerance = (std::max)(
      boost::math::tools::epsilon<RealType>(),
      (RealType)boost::math::tools::epsilon<double>() * 5) * 100;
   //
   // At float precision we need to up the tolerance, since 
   // the input values are rounded off to inexact quantities
   // the results get thrown off by a noticeable amount.
   //
   if(boost::math::tools::digits<RealType>() < 50)
      tolerance *= 50;
   if(boost::is_floating_point<RealType>::value != 1)
      tolerance *= 20; // real_concept special functions are less accurate

   cout << "Tolerance = " << tolerance << "%." << endl;

   using boost::math::chi_squared_distribution;
   using  ::boost::math::chi_squared;
   using  ::boost::math::cdf;
   using  ::boost::math::pdf;
   //
   // Test against the data from Table 6 of:
   //
   // "Self-Validating Computations of Probabilities for Selected 
   // Central and Noncentral Univariate Probability Functions."
   // Morgan C. Wang; William J. Kennedy
   // Journal of the American Statistical Association, 
   // Vol. 89, No. 427. (Sep., 1994), pp. 878-887.
   //
   test_spot(
      static_cast<RealType>(1),   // degrees of freedom
      static_cast<RealType>(6),   // non centrality
      static_cast<RealType>(0.00393),   // Chi Squared statistic
      static_cast<RealType>(0.2498463724258039e-2),       // Probability of result (CDF), P
      static_cast<RealType>(1-0.2498463724258039e-2),           // Q = 1 - P
      tolerance);
   test_spot(
      static_cast<RealType>(5),   // degrees of freedom
      static_cast<RealType>(1),   // non centrality
      static_cast<RealType>(9.23636),   // Chi Squared statistic
      static_cast<RealType>(0.8272918751175548),       // Probability of result (CDF), P
      static_cast<RealType>(1-0.8272918751175548),           // Q = 1 - P
      tolerance);
   test_spot(
      static_cast<RealType>(11),   // degrees of freedom
      static_cast<RealType>(21),   // non centrality
      static_cast<RealType>(24.72497),   // Chi Squared statistic
      static_cast<RealType>(0.2539481822183126),       // Probability of result (CDF), P
      static_cast<RealType>(1-0.2539481822183126),           // Q = 1 - P
      tolerance);
   test_spot(
      static_cast<RealType>(31),   // degrees of freedom
      static_cast<RealType>(6),   // non centrality
      static_cast<RealType>(44.98534),   // Chi Squared statistic
      static_cast<RealType>(0.8125198785064969),       // Probability of result (CDF), P
      static_cast<RealType>(1-0.8125198785064969),           // Q = 1 - P
      tolerance);
   test_spot(
      static_cast<RealType>(51),   // degrees of freedom
      static_cast<RealType>(1),   // non centrality
      static_cast<RealType>(38.56038),   // Chi Squared statistic
      static_cast<RealType>(0.8519497361859118e-1),       // Probability of result (CDF), P
      static_cast<RealType>(1-0.8519497361859118e-1),           // Q = 1 - P
      tolerance * 2);
   test_spot(
      static_cast<RealType>(100),   // degrees of freedom
      static_cast<RealType>(16),   // non centrality
      static_cast<RealType>(82.35814),   // Chi Squared statistic
      static_cast<RealType>(0.1184348822747824e-1),       // Probability of result (CDF), P
      static_cast<RealType>(1-0.1184348822747824e-1),           // Q = 1 - P
      tolerance);
   test_spot(
      static_cast<RealType>(300),   // degrees of freedom
      static_cast<RealType>(16),   // non centrality
      static_cast<RealType>(331.78852),   // Chi Squared statistic
      static_cast<RealType>(0.7355956710306709),       // Probability of result (CDF), P
      static_cast<RealType>(1-0.7355956710306709),           // Q = 1 - P
      tolerance);
   test_spot(
      static_cast<RealType>(500),   // degrees of freedom
      static_cast<RealType>(21),   // non centrality
      static_cast<RealType>(459.92612),   // Chi Squared statistic
      static_cast<RealType>(0.2797023600800060e-1),       // Probability of result (CDF), P
      static_cast<RealType>(1-0.2797023600800060e-1),           // Q = 1 - P
      tolerance);