[ceph.git] / ceph / src / boost / libs / math / test / test_nc_beta.cpp

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

//
// This must appear *before* any #includes, and precludes pch usage:
//
#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false

#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_beta.hpp> // for chi_squared_distribution
#include <boost/math/distributions/poisson.hpp> // for poisson_distribution
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp> // for test_main
#include <boost/test/results_collector.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE

#include "functor.hpp"
#include "handle_test_result.hpp"
#include "test_ncbeta_hooks.hpp"
#include "table_type.hpp"
#include "test_nc_beta.hpp"

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

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

#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
   if(boost::math::tools::digits<long double>() == 64)
   {
      //
      // Allow a small amount of error leakage from long double to double:
      //
      add_expected_result(
         "[^|]*",                          // compiler
         "[^|]*",                          // stdlib
         "[^|]*",                          // platform
         "double",                         // test type(s)
         "[^|]*large[^|]*",                // test data group
         "[^|]*", 5, 5);                   // test function
   }

   if(boost::math::tools::digits<long double>() == 64)
   {
      add_expected_result(
         "[^|]*",                          // compiler
         "[^|]*",                          // stdlib
         "[^|]*",                          // platform
         largest_type,                     // test type(s)
         "[^|]*medium[^|]*",               // test data group
         "[^|]*", 1200, 500);               // test function
      add_expected_result(
         "[^|]*",                          // compiler
         "[^|]*",                          // stdlib
         "[^|]*",                          // platform
         largest_type,                     // test type(s)
         "[^|]*large[^|]*",                // test data group
         "[^|]*", 40000, 6000);            // test function
   }
#endif
   //
   // Catch all cases come last:
   //
   add_expected_result(
      "[^|]*",                          // compiler
      "[^|]*",                          // stdlib
      "[^|]*",                          // platform
      largest_type,                     // test type(s)
      "[^|]*medium[^|]*",               // test data group
      "[^|]*", 1500, 500);               // test function
   add_expected_result(
      "[^|]*",                          // compiler
      "[^|]*",                          // stdlib
      "[^|]*",                          // platform
      "real_concept",                   // test type(s)
      "[^|]*large[^|]*",                // test data group
      "[^|]*", 30000, 4000);             // test function
   add_expected_result(
      "[^|]*",                          // compiler
      "[^|]*",                          // stdlib
      "[^|]*",                          // platform
      largest_type,                     // test type(s)
      "[^|]*large[^|]*",                // test data group
      "[^|]*", 20000, 2000);             // 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 a, RealType b, RealType lam, RealType x)
{
   using namespace boost::math;

   RealType term = pdf(poisson_distribution<RealType>(lam/2), 0)
      * ibeta_derivative(a, b, x);
   RealType sum = term;

   int i = 1;
   while(term / sum > tools::epsilon<RealType>())
   {
      term = pdf(poisson_distribution<RealType>(lam/2), i)
      * ibeta_derivative(a + i, b, x);
      ++i;
      sum += term;
   }
   return sum;
}

template <class RealType>
void test_spot(
     RealType a,     // alpha
     RealType b,     // beta
     RealType ncp,   // non-centrality param
     RealType cs,    // Chi Square statistic
     RealType P,     // CDF
     RealType Q,     // Complement of CDF
     RealType D,     // PDF
     RealType tol)   // Test tolerance
{
   boost::math::non_central_beta_distribution<RealType> dist(a, b, ncp);
   BOOST_CHECK_CLOSE(
      cdf(dist, cs), P, tol);
   //
   // Sanity checking using the naive PDF calculation above fails at
   // float precision:
   //
   if(!boost::is_same<float, RealType>::value)
   {
      BOOST_CHECK_CLOSE(
         pdf(dist, cs), naive_pdf(dist.alpha(), dist.beta(), ncp, cs), tol);
   }
   BOOST_CHECK_CLOSE(
      pdf(dist, cs), D, tol);

   if((P < 0.99) && (Q < 0.99))
   {
      //
      // We can only check this if P is not too close to 1,
      // so that we can guarantee 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);
   }
}

template <class RealType> // Any floating-point type RealType.
void test_spots(RealType)
{
   RealType tolerance = (std::max)(
      boost::math::tools::epsilon<RealType>() * 100,
      (RealType)1e-6) * 100;
   RealType abs_tolerance = boost::math::tools::epsilon<RealType>() * 100;

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

   //
   // Spot tests use values computed by the R statistical
   // package and the pbeta and dbeta functions:
   //
   test_spot(
     RealType(2),                   // alpha
     RealType(5),                   // beta
     RealType(1),                   // non-centrality param
     RealType(0.25),                // Chi Square statistic
     RealType(0.3658349),           // CDF
     RealType(1-0.3658349),         // Complement of CDF
     RealType(2.184465),            // PDF
     RealType(tolerance));
   test_spot(
     RealType(20),                  // alpha
     RealType(15),                  // beta
     RealType(35),                  // non-centrality param
     RealType(0.75),                // Chi Square statistic
     RealType(0.6994175),           // CDF
     RealType(1-0.6994175),         // Complement of CDF
     RealType(5.576146),            // PDF
     RealType(tolerance));
   test_spot(
     RealType(100),                 // alpha
     RealType(3),                   // beta
     RealType(63),                  // non-centrality param
     RealType(0.95),                // Chi Square statistic
     RealType(0.03529306),          // CDF
     RealType(1-0.03529306),        // Complement of CDF
     RealType(3.637894),            // PDF
     RealType(tolerance));
   test_spot(
     RealType(0.25),                // alpha
     RealType(0.75),                // beta
     RealType(150),                 // non-centrality param
     RealType(0.975),               // Chi Square statistic
     RealType(0.09752216),          // CDF
     RealType(1-0.09752216),        // Complement of CDF
     RealType(8.020935),            // PDF
     RealType(tolerance));

   BOOST_MATH_STD_USING
   boost::math::non_central_beta_distribution<RealType> dist(100, 3, 63);
   BOOST_CHECK_CLOSE(mean(dist), RealType(4.82280451915522329944315287538684030781836554279474240490936e13L) * exp(-RealType(31.5)) * 100 / 103, tolerance);
   // Variance only guarantees small absolute error:
   BOOST_CHECK_SMALL(variance(dist) 
      - static_cast<RealType>(RealType(4.85592267707818899235900237275021938334418424134218087127572e13L)
      * exp(RealType(-31.5)) * 100 * 101 / (103 * 104) - 
      RealType(4.82280451915522329944315287538684030781836554279474240490936e13L) * RealType(4.82280451915522329944315287538684030781836554279474240490936e13L) 
      * exp(RealType(-63)) * 10000 / (103 * 103)), abs_tolerance);
   BOOST_MATH_CHECK_THROW(skewness(dist), boost::math::evaluation_error);
   BOOST_MATH_CHECK_THROW(kurtosis(dist), boost::math::evaluation_error);
   BOOST_MATH_CHECK_THROW(kurtosis_excess(dist), boost::math::evaluation_error);
} // template <class RealType>void test_spots(RealType)


BOOST_AUTO_TEST_CASE( test_main )
{
   BOOST_MATH_CONTROL_FP;
   // Basic sanity-check spot values.
    expected_results();
   // (Parameter value, arbitrarily zero, only communicates the floating point type).
#ifdef TEST_FLOAT
   test_spots(0.0F); // Test float.
#endif
#ifdef TEST_DOUBLE
   test_spots(0.0); // Test double.
#endif
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
#ifdef TEST_LDOUBLE
   test_spots(0.0L); // Test long double.
#endif
#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
#ifdef TEST_REAL_CONCEPT
   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
#endif
#endif
#endif

#ifdef TEST_FLOAT
   test_accuracy(0.0F, "float"); // Test float.
#endif
#ifdef TEST_DOUBLE
   test_accuracy(0.0, "double"); // Test double.
#endif
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
#ifdef TEST_LDOUBLE
   test_accuracy(0.0L, "long double"); // Test long double.
#endif
#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
#ifdef TEST_REAL_CONCEPT
   test_accuracy(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
#endif
#endif
#endif
   
} // BOOST_AUTO_TEST_CASE( test_main )
Commit	Line	Data
7c673cae FG	1	// test_nc_beta.cpp
	2
	3	// Copyright John Maddock 2008.
	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	//
	11	// This must appear before any #includes, and precludes pch usage:
	12	//
	13	#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
	14
	15	#ifdef _MSC_VER
	16	#pragma warning (disable:4127 4512)
	17	#endif
	18
	19	#if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
	20	# define TEST_FLOAT
	21	# define TEST_DOUBLE
	22	# define TEST_LDOUBLE
	23	# define TEST_REAL_CONCEPT
	24	#endif
	25
	26	#include <boost/math/concepts/real_concept.hpp> // for real_concept
	27	#include <boost/math/distributions/non_central_beta.hpp> // for chi_squared_distribution
	28	#include <boost/math/distributions/poisson.hpp> // for poisson_distribution
	29	#define BOOST_TEST_MAIN
	30	#include <boost/test/unit_test.hpp> // for test_main
	31	#include <boost/test/results_collector.hpp>
	32	#include <boost/test/unit_test.hpp>
92f5a8d4	33	#include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
7c673cae FG	34
	35	#include "functor.hpp"
	36	#include "handle_test_result.hpp"
	37	#include "test_ncbeta_hooks.hpp"
	38	#include "table_type.hpp"
	39	#include "test_nc_beta.hpp"
	40
	41	#include <iostream>
	42	using std::cout;
	43	using std::endl;
	44	#include <limits>
	45	using std::numeric_limits;
	46
	47	void expected_results()
	48	{
	49	//
	50	// Define the max and mean errors expected for
	51	// various compilers and platforms.
	52	//
	53	const char* largest_type;
	54	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
	55	if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
	56	{
	57	largest_type = "(long\\s+)?double\|real_concept";
	58	}
	59	else
	60	{
	61	largest_type = "long double\|real_concept";
	62	}
	63	#else
	64	largest_type = "(long\\s+)?double\|real_concept";
	65	#endif
	66
	67	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
	68	if(boost::math::tools::digits<long double>() == 64)
	69	{
	70	//
	71	// Allow a small amount of error leakage from long double to double:
	72	//
	73	add_expected_result(
	74	"[^\|]*", // compiler
	75	"[^\|]*", // stdlib
	76	"[^\|]*", // platform
	77	"double", // test type(s)
	78	"[^\|]large[^\|]", // test data group
	79	"[^\|]*", 5, 5); // test function
	80	}
	81
	82	if(boost::math::tools::digits<long double>() == 64)
	83	{
	84	add_expected_result(
	85	"[^\|]*", // compiler
	86	"[^\|]*", // stdlib
	87	"[^\|]*", // platform
	88	largest_type, // test type(s)
	89	"[^\|]medium[^\|]", // test data group
	90	"[^\|]*", 1200, 500); // test function
	91	add_expected_result(
	92	"[^\|]*", // compiler
	93	"[^\|]*", // stdlib
	94	"[^\|]*", // platform
	95	largest_type, // test type(s)
	96	"[^\|]large[^\|]", // test data group
	97	"[^\|]*", 40000, 6000); // test function
98	}
99	#endif
100	//
101	// Catch all cases come last:
102	//
103	add_expected_result(
104	"[^\|]*", // compiler
105	"[^\|]*", // stdlib
106	"[^\|]*", // platform
107	largest_type, // test type(s)
108	"[^\|]medium[^\|]", // test data group
109	"[^\|]*", 1500, 500); // test function
110	add_expected_result(
111	"[^\|]*", // compiler
112	"[^\|]*", // stdlib
113	"[^\|]*", // platform
114	"real_concept", // test type(s)
115	"[^\|]large[^\|]", // test data group
116	"[^\|]*", 30000, 4000); // test function
117	add_expected_result(
118	"[^\|]*", // compiler
119	"[^\|]*", // stdlib
120	"[^\|]*", // platform
121	largest_type, // test type(s)
122	"[^\|]large[^\|]", // test data group
123	"[^\|]*", 20000, 2000); // test function
124	//
125	// Finish off by printing out the compiler/stdlib/platform names,
126	// we do this to make it easier to mark up expected error rates.
127	//
128	std::cout << "Tests run with " << BOOST_COMPILER << ", "
129	<< BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
130	}
131
132	template <class RealType>
133	RealType naive_pdf(RealType a, RealType b, RealType lam, RealType x)
134	{
135	using namespace boost::math;
136
137	RealType term = pdf(poisson_distribution<RealType>(lam/2), 0)
138	* ibeta_derivative(a, b, x);
139	RealType sum = term;
140
141	int i = 1;
142	while(term / sum > tools::epsilon<RealType>())
143	{
144	term = pdf(poisson_distribution<RealType>(lam/2), i)
145	* ibeta_derivative(a + i, b, x);
146	++i;
147	sum += term;
148	}
149	return sum;
150	}
151
152	template <class RealType>
153	void test_spot(
154	RealType a, // alpha
155	RealType b, // beta
156	RealType ncp, // non-centrality param
157	RealType cs, // Chi Square statistic
158	RealType P, // CDF
159	RealType Q, // Complement of CDF
160	RealType D, // PDF
161	RealType tol) // Test tolerance
162	{
163	boost::math::non_central_beta_distribution<RealType> dist(a, b, ncp);
164	BOOST_CHECK_CLOSE(
165	cdf(dist, cs), P, tol);
166	//
167	// Sanity checking using the naive PDF calculation above fails at
168	// float precision:
169	//
170	if(!boost::is_same<float, RealType>::value)
171	{
172	BOOST_CHECK_CLOSE(
173	pdf(dist, cs), naive_pdf(dist.alpha(), dist.beta(), ncp, cs), tol);
174	}
175	BOOST_CHECK_CLOSE(
176	pdf(dist, cs), D, tol);
177
178	if((P < 0.99) && (Q < 0.99))
179	{
180	//
181	// We can only check this if P is not too close to 1,
182	// so that we can guarantee Q is reasonably free of error:
183	//
184	BOOST_CHECK_CLOSE(
185	cdf(complement(dist, cs)), Q, tol);
186	BOOST_CHECK_CLOSE(
187	quantile(dist, P), cs, tol * 10);
188	BOOST_CHECK_CLOSE(
189	quantile(complement(dist, Q)), cs, tol * 10);
190	}
191	}
192
193	template <class RealType> // Any floating-point type RealType.
194	void test_spots(RealType)
195	{
196	RealType tolerance = (std::max)(
197	boost::math::tools::epsilon<RealType>() * 100,
198	(RealType)1e-6) * 100;
199	RealType abs_tolerance = boost::math::tools::epsilon<RealType>() * 100;
200
201	cout << "Tolerance = " << tolerance << "%." << endl;
202
203	//
204	// Spot tests use values computed by the R statistical
205	// package and the pbeta and dbeta functions:
206	//
207	test_spot(
208	RealType(2), // alpha
209	RealType(5), // beta
210	RealType(1), // non-centrality param
211	RealType(0.25), // Chi Square statistic
212	RealType(0.3658349), // CDF
213	RealType(1-0.3658349), // Complement of CDF
214	RealType(2.184465), // PDF
215	RealType(tolerance));
216	test_spot(
217	RealType(20), // alpha
218	RealType(15), // beta
219	RealType(35), // non-centrality param
220	RealType(0.75), // Chi Square statistic
221	RealType(0.6994175), // CDF
222	RealType(1-0.6994175), // Complement of CDF
223	RealType(5.576146), // PDF
224	RealType(tolerance));
225	test_spot(
226	RealType(100), // alpha
227	RealType(3), // beta
228	RealType(63), // non-centrality param
229	RealType(0.95), // Chi Square statistic
230	RealType(0.03529306), // CDF
231	RealType(1-0.03529306), // Complement of CDF
232	RealType(3.637894), // PDF
233	RealType(tolerance));
234	test_spot(
235	RealType(0.25), // alpha
236	RealType(0.75), // beta
237	RealType(150), // non-centrality param
238	RealType(0.975), // Chi Square statistic
239	RealType(0.09752216), // CDF
240	RealType(1-0.09752216), // Complement of CDF
241	RealType(8.020935), // PDF
242	RealType(tolerance));
243
244	BOOST_MATH_STD_USING
245	boost::math::non_central_beta_distribution<RealType> dist(100, 3, 63);
246	BOOST_CHECK_CLOSE(mean(dist), RealType(4.82280451915522329944315287538684030781836554279474240490936e13L) * exp(-RealType(31.5)) * 100 / 103, tolerance);
247	// Variance only guarantees small absolute error:
248	BOOST_CHECK_SMALL(variance(dist)
249	- static_cast<RealType>(RealType(4.85592267707818899235900237275021938334418424134218087127572e13L)
250	* exp(RealType(-31.5)) * 100 * 101 / (103 * 104) -
251	RealType(4.82280451915522329944315287538684030781836554279474240490936e13L) * RealType(4.82280451915522329944315287538684030781836554279474240490936e13L)
252	* exp(RealType(-63)) * 10000 / (103 * 103)), abs_tolerance);
253	BOOST_MATH_CHECK_THROW(skewness(dist), boost::math::evaluation_error);
254	BOOST_MATH_CHECK_THROW(kurtosis(dist), boost::math::evaluation_error);
255	BOOST_MATH_CHECK_THROW(kurtosis_excess(dist), boost::math::evaluation_error);
256	} // template <class RealType>void test_spots(RealType)
257
258
259	BOOST_AUTO_TEST_CASE( test_main )
260	{
261	BOOST_MATH_CONTROL_FP;
262	// Basic sanity-check spot values.
263	expected_results();
264	// (Parameter value, arbitrarily zero, only communicates the floating point type).
265	#ifdef TEST_FLOAT
266	test_spots(0.0F); // Test float.
267	#endif
268	#ifdef TEST_DOUBLE
269	test_spots(0.0); // Test double.
270	#endif
271	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
272	#ifdef TEST_LDOUBLE
273	test_spots(0.0L); // Test long double.
274	#endif
275	#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
276	#ifdef TEST_REAL_CONCEPT
277	test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
278	#endif
279	#endif
280	#endif
281
282	#ifdef TEST_FLOAT
283	test_accuracy(0.0F, "float"); // Test float.
284	#endif
285	#ifdef TEST_DOUBLE
286	test_accuracy(0.0, "double"); // Test double.
287	#endif
288	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
289	#ifdef TEST_LDOUBLE
290	test_accuracy(0.0L, "long double"); // Test long double.
291	#endif
292	#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
293	#ifdef TEST_REAL_CONCEPT
294	test_accuracy(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
295	#endif
296	#endif
297	#endif
298
299	} // BOOST_AUTO_TEST_CASE( test_main )
300