boost/random
2
0
Fork 0
mirror of https://github.com/boostorg/random.git synced 2026-01-19 04:22:17 +00:00
Code Issues Projects Releases Wiki Activity

This commit was manufactured by cvs2svn to create branch 'conversion'.

Browse Source 
[SVN r7663]
This commit is contained in:
nobody
2000-07-30 10:33:54 +00:00
parent 28cad88d43
commit a130d643b4
9 changed files with 0 additions and 2301 deletions
Download Patch File Download Diff File Expand all files Collapse all files

147
histogram.cpp
View File

@@ -1,147 +0,0 @@
/* boost histogram.cpp graphical verification of distribution functions
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
* This test program allows to visibly examine the results of the
* distribution functions.
*/
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <cmath>
#include <string>
#include <boost/random.hpp>
void plot_histogram(const std::vector<int>& slots, int samples,
double from, double to)
{
int m = *std::max_element(slots.begin(), slots.end());
const int nRows = 20;
std::cout.setf(std::ios::fixed|std::ios::left);
std::cout.precision(5);
for(int r = 0; r < nRows; r++) {
double y = ((nRows - r) * double(m))/(nRows * samples);
std::cout << std::setw(10) << y << " ";
for(unsigned int col = 0; col < slots.size(); col++) {
char out = ' ';
if(slots[col]/double(samples) >= y)
out = 'x';
std::cout << out;
}
std::cout << std::endl;
}
std::cout << std::setw(12) << " "
<< std::setw(10) << from;
std::cout.setf(std::ios::right, std::ios::adjustfield);
std::cout << std::setw(slots.size()-10) << to << std::endl;
}
// I am not sure whether these two should be in the library as well
// maintain sum of NumberGenerator results
template<class NumberGenerator,
class Sum = typename NumberGenerator::result_type>
class sum_result
{
public:
typedef NumberGenerator base_type;
typedef typename base_type::result_type result_type;
explicit sum_result(const base_type & g) : gen(g), _sum(0) { }
result_type operator()() { result_type r = gen(); _sum += r; return r; }
base_type & base() { return gen; }
Sum sum() const { return _sum; }
void reset() { _sum = 0; }
private:
base_type gen;
Sum _sum;
};
// maintain square sum of NumberGenerator results
template<class NumberGenerator,
class Sum = typename NumberGenerator::result_type>
class squaresum_result
{
public:
typedef NumberGenerator base_type;
typedef typename base_type::result_type result_type;
explicit squaresum_result(const base_type & g) : gen(g), _sum(0) { }
result_type operator()() { result_type r = gen(); _sum += r*r; return r; }
base_type & base() { return gen; }
Sum squaresum() const { return _sum; }
void reset() { _sum = 0; }
private:
base_type gen;
Sum _sum;
};
template<class RNG>
void histogram(RNG base, int samples, double from, double to,
const std::string & name)
{
typedef squaresum_result<sum_result<RNG, double>, double > SRNG;
SRNG gen((sum_result<RNG, double>(base)));
const int nSlots = 60;
std::vector<int> slots(nSlots,0);
for(int i = 0; i < samples; i++) {
double val = gen();
if(val < from || val >= to) // early check avoids overflow
continue;
int slot = int((val-from)/(to-from) * nSlots);
if(slot < 0 || slot > (int)slots.size())
continue;
slots[slot]++;
}
std::cout << name << std::endl;
plot_histogram(slots, samples, from, to);
double mean = gen.base().sum() / samples;
std::cout << "mean: " << mean
<< " sigma: " << std::sqrt(gen.squaresum()/samples-mean*mean)
<< "\n" << std::endl;
}
template<class PRNG>
void histograms()
{
PRNG rng;
using namespace boost;
histogram(uniform_01<PRNG>(rng), 100000, -0.5, 1.5, "uniform_01(0,1)");
histogram(bernoulli_distribution<PRNG>(rng, 0.2), 100000, -0.5, 1.5,
"bernoulli(0.2)");
histogram(triangle_distribution<PRNG>(rng, 1, 2, 8), 100000, 0, 10,
"triangle(1,2,8)");
histogram(geometric_distribution<PRNG>(rng, 5.0/6.0), 100000, 0, 10,
"geometric(5/6)");
histogram(exponential_distribution<PRNG>(rng, 0.3), 100000, 0, 10,
"exponential(0.3)");
histogram(cauchy_distribution<PRNG>(rng), 100000, -5, 5,
"cauchy");
histogram(lognormal_distribution<PRNG>(rng, 3, 2), 100000, 0, 10,
"lognormal");
histogram(normal_distribution<PRNG>(rng), 100000, -3, 3,
"normal");
histogram(normal_distribution<PRNG>(rng, 0.5, 0.5), 100000, -3, 3,
"normal(0.5, 0.5)");
}
int main()
{
histograms<boost::mt19937>();
}

83
integrate.hpp
View File

@@ -1,83 +0,0 @@
/* integrate.hpp header file
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
* Revision history
*/
#ifndef INTEGRATE_HPP
#define INTEGRATE_HPP
#include <limits>
template<class UnaryFunction>
inline typename UnaryFunction::result_type
trapezoid(UnaryFunction f, typename UnaryFunction::argument_type a,
typename UnaryFunction::argument_type b, int n)
{
typename UnaryFunction::result_type tmp = 0;
for(int i = 1; i <= n-1; ++i)
tmp += f(a+(b-a)/n*i);
return (b-a)/2/n * (f(a) + f(b) + 2*tmp);
}
template<class UnaryFunction>
inline typename UnaryFunction::result_type
simpson(UnaryFunction f, typename UnaryFunction::argument_type a,
typename UnaryFunction::argument_type b, int n)
{
typename UnaryFunction::result_type tmp1 = 0;
for(int i = 1; i <= n-1; ++i)
tmp1 += f(a+(b-a)/n*i);
typename UnaryFunction::result_type tmp2 = 0;
for(int i = 1; i <= n ; ++i)
tmp2 += f(a+(b-a)/2/n*(2*i-1));
return (b-a)/6/n * (f(a) + f(b) + 2*tmp1 + 4*tmp2);
}
// compute b so that f(b) = y; assume f is monotone increasing
template<class UnaryFunction>
inline typename UnaryFunction::argument_type
invert_monotone_inc(UnaryFunction f, typename UnaryFunction::result_type y,
typename UnaryFunction::argument_type lower = -1,
typename UnaryFunction::argument_type upper = 1)
{
while(upper-lower > 1e-6) {
double middle = (upper+lower)/2;
if(f(middle) > y)
upper = middle;
else
lower = middle;
}
return (upper+lower)/2;
}
// compute b so that I(f(x), a, b) == y
template<class UnaryFunction>
inline typename UnaryFunction::argument_type
quantil(UnaryFunction f, typename UnaryFunction::argument_type a,
typename UnaryFunction::result_type y,
typename UnaryFunction::argument_type step)
{
typedef typename UnaryFunction::result_type result_type;
if(y >= 1.0)
return std::numeric_limits<result_type>::infinity();
typename UnaryFunction::argument_type b = a;
for(result_type result = 0; result < y; b += step)
result += step*f(b);
return b;
}
#endif /* INTEGRATE_HPP */

67
nondet_random_speed.cpp
View File

@@ -1,67 +0,0 @@
/* boost nondet_random_speed.cpp performance test
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
*/
#include <iostream>
#include <string>
#include <boost/timer.hpp>
#include <boost/nondet_random.hpp>
// set to your CPU frequency in MHz
static const double cpu_frequency = 200 * 1e6;
static void show_elapsed(double end, int iter, const std::string & name)
{
double usec = end/iter*1e6;
double cycles = usec * cpu_frequency/1e6;
std::cout << name << ": "
<< usec*1e3 << " nsec/loop = "
<< cycles << " CPU cycles"
<< std::endl;
}
template<class Result, class RNG>
static void timing(RNG & rng, int iter, const std::string& name)
{
volatile Result tmp; // make sure we're not optimizing too much
boost::timer t;
for(int i = 0; i < iter; i++)
tmp = rng();
show_elapsed(t.elapsed(), iter, name);
}
template<class RNG>
void run(int iter, const std::string & name)
{
RNG rng;
timing<long>(rng, iter, name);
}
int main(int argc, char*argv[])
{
if(argc != 2) {
std::cerr << "usage: " << argv[0] << " iterations" << std::endl;
return 1;
}
int iter = std::atoi(argv[1]);
#ifdef __linux__
boost::random_device dev;
timing<unsigned int>(dev, iter, "random_device");
#endif
return 0;
}

80
random_demo.cpp
View File

@@ -1,80 +0,0 @@
/* boost random_demo.cpp profane demo
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
* A short demo program how to use the random number library.
*/
#include <iostream>
#include <fstream>
#include <ctime> // std::time
#include <boost/random.hpp>
// try boost::mt19937 or boost::ecuyer1988 instead of boost::minstd_rand
typedef boost::minstd_rand base_generator_type;
// This is a reproducible simulation experiment.
void experiment(base_generator_type & generator)
{
boost::uniform_smallint<base_generator_type> die(generator, 1, 6);
// you can use a STL Iterator interface
for(int i = 0; i < 10; i++)
std::cout << *die++ << " ";
std::cout << '\n';
}
int main()
{
// initialize by reproducible seed
base_generator_type generator(42);
std::cout << "10 samples of a uniform distribution in [0..1):\n";
boost::uniform_01<base_generator_type> uni(generator);
std::cout.setf(std::ios::fixed);
// you can also use a STL Generator interface
for(int i = 0; i < 10; i++)
std::cout << uni() << '\n';
// change seed to something else
// Note: this is not the preferred way of hacking around missing std::
generator.seed(
#ifndef BOOST_NO_STDC_NAMESPACE
std::
#endif
time(0));
std::cout << "\nexperiment: roll a die 10 times:\n";
base_generator_type saved_generator = generator;
experiment(generator);
std::cout << "redo the experiment to verify it:\n";
experiment(saved_generator);
// after that, both generators are equivalent
assert(generator == saved_generator);
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
{
// save the generator state for future use,
// can be read again at any time via operator>>
std::ofstream file("rng.saved", std::ofstream::trunc);
file << generator;
}
#endif
// Some compilers don't pay attention to std:3.6.1/5 and issue a
// warning here if "return 0;" is omitted.
return 0;
}

108
random_device.cpp
View File

@@ -1,108 +0,0 @@
/* boost random_device.cpp implementation
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
*/
#include <boost/nondet_random.hpp>
#include <string>
#ifdef __linux__
// the default is the unlimited capacity device, using some secure hash
// try "/dev/random" for blocking when the entropy pool has drained
const char * const boost::random_device::default_token = "/dev/urandom";
/*
* This uses the POSIX interface for unbuffered reading.
* Using buffered std::istream would consume entropy which may
* not actually be used. Entropy is a precious good we avoid
* wasting.
*/
#if defined(__GNUC__) && defined(_CXXRT_STD_NAME)
// I have severe difficulty to get the POSIX includes to work with
// -fhonor-std and Dietmar Kühl's standard C++ library. Hack around that
// problem for now.
extern "C" {
static const int O_RDONLY = 0;
extern int open(const char *__file, int __oflag, ...);
extern int read(int __fd, __ptr_t __buf, size_t __nbytes);
extern int close(int __fd);
}
#else
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h> // open
#include <unistd.h> // read, close
#endif
#include <errno.h> // errno
#include <string.h> // strerror
#include <stdexcept> // std::invalid_argument
class boost::random_device::impl
{
public:
impl(const std::string & token) : path(token) {
fd = open(token.c_str(), O_RDONLY);
if(fd < 0)
error("cannot open");
}
~impl() { if(close(fd) < 0) error("could not close"); }
unsigned int next() {
unsigned int result;
long sz = read(fd, reinterpret_cast<char *>(&result), sizeof(result));
if(sz == -1)
error("error while reading");
else if(sz != sizeof(result)) {
errno = 0;
error("EOF while reading");
}
return result;
}
private:
void error(const std::string & msg) {
throw std::invalid_argument("boost::random_device: " + msg +
" random-number pseudo-device " + path +
": " + strerror(errno));
}
const std::string path;
int fd;
};
#endif // __linux__
boost::random_device::random_device(const std::string& token)
: pimpl(new impl(token))
{
assert(std::numeric_limits<result_type>::max() == max_value);
}
boost::random_device::~random_device()
{
// the complete class impl is now visible, so we're safe
// (see comment in random.hpp)
delete pimpl;
}
unsigned int boost::random_device::operator()()
{
return pimpl->next();
}

217
random_speed.cpp
View File

@@ -1,217 +0,0 @@
/* boost random_speed.cpp performance measurements
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*/
#include <iostream>
#include <cstdlib>
#include <string>
#include <boost/random.hpp>
#include <boost/progress.hpp>
/*
* Configuration Section
*/
// define if your C library supports the non-standard drand48 family
#undef HAVE_DRAND48
// define if you have the original mt19937int.c (with commented out main())
#undef HAVE_MT19937INT_C
// set to your CPU frequency in MHz
static const double cpu_frequency = 200 * 1e6;
/*
* End of Configuration Section
*/
/*
* General portability note:
* MSVC mis-compiles explicit function template instantiations.
* For example, f<A>() and f<B>() are both compiled to call f<A>().
* BCC is unable to implicitly convert a "const char *" to a std::string
* when using explicit function template instantiations.
*
* Therefore, avoid explicit function template instantiations.
*/
// provides a run-time configurable linear congruential generator, just
// for comparison
template<class IntType>
class linear_congruential
{
public:
typedef IntType result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = false;
#else
enum { has_fixed_range = false };
#endif
linear_congruential(IntType x0, IntType a, IntType c, IntType m)
: _x(x0), _a(a), _c(c), _m(m) { }
// compiler-generated copy ctor and assignment operator are fine
void seed(IntType x0, IntType a, IntType c, IntType m)
{ _x = x0; _a = a; _c = c; _m = m; }
void seed(IntType x0) { _x = x0; }
result_type operator()() { _x = (_a*_x+_c) % _m; return _x; }
result_type min() const { return _c == 0 ? 1 : 0; }
result_type max() const { return _m -1; }
private:
IntType _x, _a, _c, _m;
};
static void show_elapsed(double end, int iter, const std::string & name)
{
double usec = end/iter*1e6;
double cycles = usec * cpu_frequency/1e6;
std::cout << name << ": "
<< usec*1e3 << " nsec/loop = "
<< cycles << " CPU cycles"
<< std::endl;
}
template<class Result, class RNG>
static void timing(RNG & rng, int iter, const std::string& name, const Result&)
{
volatile Result tmp; // make sure we're not optimizing too much
boost::timer t;
for(int i = 0; i < iter; i++)
tmp = rng();
show_elapsed(t.elapsed(), iter, name);
}
template<class RNG>
static void timing_sphere(RNG & rng, int iter, const std::string & name)
{
boost::timer t;
for(int i = 0; i < iter; i++) {
// the special return value convention of uniform_on_sphere saves 20% CPU
const std::vector<double> & tmp = rng();
(void) tmp[0];
}
show_elapsed(t.elapsed(), iter, name);
}
template<class RNG>
void run(int iter, const std::string & name, const RNG &)
{
RNG rng;
std::cout << (RNG::has_fixed_range ? "fixed-range " : "");
// BCC has trouble with string autoconversion for explicit specializations
timing(rng, iter, std::string(name), 0l);
}
#ifdef HAVE_DRAND48
// requires non-standard C library support for srand48/lrand48
void run(int iter, const std::string & name, int)
{
std::srand48(1);
timing(std::lrand48, iter, name, 0l);
}
#endif
#ifdef HAVE_MT19937INT_C // requires the original mt19937int.c
extern "C" void sgenrand(unsigned long);
extern "C" unsigned long genrand();
void run(int iter, const std::string & name, float)
{
sgenrand(4357);
timing(genrand, iter, name, 0u);
}
#endif
template<class Gen>
void distrib(int iter, const std::string & name, const Gen &)
{
Gen gen;
boost::uniform_smallint<Gen> usmallint(gen, -2, 4);
timing(usmallint, iter, name + " uniform_smallint", 0);
boost::uniform_int<Gen> uint(gen, -2, 4);
timing(uint, iter, name + " uniform_int", 0);
boost::uniform_01<Gen> uni(gen);
timing(uni, iter, name + " uniform_01", 0.0);
boost::uniform_real<Gen> ur(gen, -5.3, 4.8);
timing(ur, iter, name + " uniform_real", 0.0);
boost::geometric_distribution<Gen> geo(gen, 0.5);
timing(geo, iter, name + " geometric", 0);
boost::triangle_distribution<Gen> tria(gen, 1, 2, 7);
timing(tria, iter, name + " triangle", 0.0);
boost::exponential_distribution<Gen> ex(gen, 3);
timing(ex, iter, name + " exponential", 0.0);
boost::normal_distribution<Gen,double> no2(gen);
timing(no2, iter, name + " normal polar", 0.0);
boost::lognormal_distribution<Gen,double> lnorm(gen, 1, 1);
timing(lnorm, iter, name + " lognormal", 0.0);
boost::uniform_on_sphere<Gen> usph(gen, 3);
timing_sphere(usph, iter/10, name + " uniform_on_sphere");
}
int main(int argc, char*argv[])
{
if(argc != 2) {
std::cerr << "usage: " << argv[0] << " iterations" << std::endl;
return 1;
}
// okay, it's ugly, but it's only used here
int iter =
#ifndef BOOST_NO_STDC_NAMESPACE
std::
#endif
atoi(argv[1]);
#ifdef BOOST_STDINT_H_HAS_UINT64_T
run(iter, "rand48", boost::rand48);
linear_congruential<unsigned long long> lcg48(1ULL<<16 | 0x330e,
0x5DEECE66DULL, 0xB,
(1ULL<<48));
timing(lcg48, iter, "lrand48 run-time", 0l);
#endif
#ifdef HAVE_DRAND48
// requires non-standard C library support for srand48/lrand48
run(iter, "lrand48", 0); // coded for lrand48()
#endif
run(iter, "minstd_rand", boost::minstd_rand());
run(iter, "ecuyer combined", boost::ecuyer1988());
run(iter, "kreutzer1986", boost::kreutzer1986());
run(iter, "hellekalek1995 (inversive)", boost::hellekalek1995());
run(iter, "mt11213b", boost::mt11213b());
run(iter, "mt19937", boost::mt19937());
#ifdef HAVE_MT19937INT_C
// requires the original mt19937int.c
run<float>(iter, "mt19937 original"); // coded for sgenrand()/genrand()
#endif
distrib(iter, "minstd_rand", boost::minstd_rand());
distrib(iter, "kreutzer1986", boost::kreutzer1986());
distrib(iter, "mt19937", boost::mt19937());
};

309
random_test.cpp
View File

@@ -1,309 +0,0 @@
/* boost random_test.cpp various tests
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*/
#include <iostream>
#include <fstream>
#include <string>
#include <cassert>
#include <boost/random.hpp>
/*
* General portability note:
* MSVC mis-compiles explicit function template instantiations.
* For example, f<A>() and f<B>() are both compiled to call f<A>().
* BCC is unable to implicitly convert a "const char *" to a std::string
* when using explicit function template instantiations.
*
* Therefore, avoid explicit function template instantiations.
*/
/*
* Validate correct implementation
*/
template<class PRNG>
void validate(const std::string & name, const PRNG &)
{
std::cout << "validating " << name << ": ";
PRNG rng;
for(int i = 0; i < 9999; i++)
rng();
typename PRNG::result_type val = rng();
bool result = const_cast<const PRNG&>(rng).validation(val);
// allow for a simple eyeball check for MSVC instantiation brokenness
// (if the numbers for all generators are the same, it's obviously broken)
std::cout << val << " ";
if(!result)
std::cout << "validation failed or not possible.\n";
else
std::cout << "ok\n";
}
void validate_all()
{
using namespace boost;
#ifdef BOOST_STDINT_H_HAS_UINT64_T
validate("rand48", rand48());
#endif
validate("minstd_rand", minstd_rand());
validate("minstd_rand0", minstd_rand0());
validate("ecuyer combined", ecuyer1988());
validate("mt19937", mt19937());
validate("kreutzer1986", kreutzer1986());
}
/*
* Check function signatures
*/
template<class Generator>
void instantiate_iterator_interface(Generator & gen)
{
++gen;
gen++;
typename Generator::result_type res = *gen;
typename Generator::pointer p = &res;
(void) p;
typename Generator::reference ref = res;
(void) ref;
typename Generator::difference_type diff = 0;
(void) diff;
std::input_iterator_tag it(typename Generator::iterator_category());
(void) it;
assert(res == *gen++);
assert(gen == gen);
Generator gen2 = gen;
assert(gen == gen2); // must be equal to a copy
#if 0
++gen2;
// this is only correct for elementary generators, others have ref members
assert(gen != gen2);
#endif
}
template<class URNG, class ResultType>
void instantiate_urng(const std::string & s, const URNG &, const ResultType &)
{
std::cout << "Basic tests for " << s << std::endl;
URNG urng;
int a[URNG::has_fixed_range ? 5 : 10]; // compile-time constant
(void) a; // avoid "unused" warning
typename URNG::result_type x1 = urng();
ResultType x2 = x1;
(void) x2; // avoid "unused" warning
#ifndef BOOST_MSVC // MSVC brokenness
URNG urng2(urng); // copy constructor
assert(urng == urng2); // operator==
urng();
urng2 = urng; // assignment
assert(urng == urng2);
#endif // BOOST_MSVC
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
// Streamable concept not supported for broken compilers
{
std::ofstream file("rng.tmp", std::ofstream::trunc);
file << urng;
}
// move forward
urng();
{
// restore old state
std::ifstream file("rng.tmp");
file >> urng;
}
#ifndef BOOST_MSVC // MSVC brokenness
assert(urng == urng2);
#endif // BOOST_MSVC
#endif // BOOST_NO_OPERATORS_IN_NAMESPACE
// instantiate various distributions with this URNG
boost::uniform_smallint<URNG> unismall(urng, 0, 11);
instantiate_iterator_interface(unismall);
unismall();
boost::uniform_int<URNG> uni_int(urng, -200, 20000);
instantiate_iterator_interface(uni_int);
uni_int();
boost::uniform_real<URNG> uni_real(urng, 0, 2.1);
instantiate_iterator_interface(uni_real);
uni_real();
boost::bernoulli_distribution<URNG> ber(urng, 0.2);
instantiate_iterator_interface(ber);
ber();
boost::geometric_distribution<URNG> geo(urng, 0.8);
instantiate_iterator_interface(geo);
geo();
boost::triangle_distribution<URNG> tria(urng, 1, 1.5, 7);
instantiate_iterator_interface(tria);
tria();
boost::exponential_distribution<URNG> ex(urng, 5);
instantiate_iterator_interface(ex);
ex();
boost::normal_distribution<URNG> norm(urng);
instantiate_iterator_interface(norm);
norm();
boost::lognormal_distribution<URNG> lnorm(urng, 1, 1);
instantiate_iterator_interface(lnorm);
lnorm();
boost::uniform_on_sphere<URNG> usph(urng, 2);
instantiate_iterator_interface(usph);
usph();
}
void instantiate_all()
{
using namespace boost;
#ifdef BOOST_STDINT_H_HAS_UINT64_T
instantiate_urng("rand48", rand48, 0);
rand48 rnd(5);
rand48 rnd2(uint64_t(0x80000000) * 42);
rnd.seed(17);
rnd2.seed(uint64_t(0x80000000) * 49);
#endif
instantiate_urng("minstd_rand0", minstd_rand0(), 0);
instantiate_urng("minstd_rand", minstd_rand(), 0);
minstd_rand mstd(42);
mstd.seed(17);
instantiate_urng("ecuyer1988", ecuyer1988(), 0);
instantiate_urng("kreutzer1986", kreutzer1986(), 0);
instantiate_urng("hellekalek1995", hellekalek1995(), 0);
instantiate_urng("mt11213b", mt11213b(), 0u);
instantiate_urng("mt19937", mt19937(), 0u);
mt19937 mt(17u); // needs to be an exact type match for MSVC
mt.seed(42u);
mt19937 mt2(mstd);
mt2.seed(mstd);
random_number_generator<mt19937> std_rng(mt2);
(void) std_rng(10);
}
/*
* A few equidistribution tests
*/
// yet to come...
template<class Generator>
void check_uniform_int(Generator & gen, int iter)
{
std::cout << "testing uniform_int(" << gen.min() << "," << gen.max()
<< ")" << std::endl;
int range = gen.max()-gen.min()+1;
std::vector<int> bucket(range);
for(int j = 0; j < iter; j++) {
int result = gen();
if(result < gen.min() || result > gen.max())
std::cerr << " ... delivers " << result << std::endl;
else
bucket[result-gen.min()]++;
}
int sum = 0;
// use a different variable name "k", because MSVC has broken "for" scoping
for(int k = 0; k < range; k++)
sum += bucket[k];
int avg = sum/range;
for(int i = 0; i < range; i++) {
if(abs(bucket[i] - avg) > 2*avg/sqrt(iter)) { // 95% confidence interval
std::cout << " ... has bucket[" << i << "] = " << bucket[i]
<< " (distance " << (bucket[i] - avg) << ")"
<< std::endl;
}
}
}
template<class Generator>
void test_uniform_int(Generator & gen)
{
typedef boost::uniform_int<Generator, int> int_gen;
// large range => small range (modulo case)
int_gen uint12(gen,1,2);
assert(uint12.min() == 1);
assert(uint12.max() == 2);
check_uniform_int(uint12, 100000);
int_gen uint16(gen,1,6);
check_uniform_int(uint16, 100000);
// test chaining to get all cases in operator()
typedef boost::uniform_int<int_gen, int> intint_gen;
// identity map
intint_gen uint01(uint12, 0, 1);
check_uniform_int(uint01, 100000);
// small range => larger range
intint_gen uint05(uint12, -3, 2);
check_uniform_int(uint05, 100000);
typedef boost::uniform_int<intint_gen, int> intintint_gen;
// small => larger range, not power of two
// (for unknown reasons, this has noticeably uneven distribution)
intintint_gen uint1_49(uint05, 1, 49);
check_uniform_int(uint1_49, 500000);
// larger => small range, rejection case
intintint_gen uint1_4(uint05, 1, 4);
check_uniform_int(uint1_4, 100000);
}
#if defined(BOOST_MSVC) && _MSC_VER <= 1200
// These explicit instantiations are necessary, otherwise MSVC does
// find the boost/operators.hpp inline friends.
// We ease the typing with a suitable preprocessor macro.
#define INSTANT(x) \
template class boost::uniform_smallint<x>; \
template class boost::uniform_int<x>; \
template class boost::uniform_real<x>; \
template class boost::bernoulli_distribution<x>; \
template class boost::geometric_distribution<x>; \
template class boost::triangle_distribution<x>; \
template class boost::exponential_distribution<x>; \
template class boost::normal_distribution<x>; \
template class boost::uniform_on_sphere<x>; \
template class boost::lognormal_distribution<x>;
INSTANT(boost::minstd_rand0)
INSTANT(boost::minstd_rand)
INSTANT(boost::ecuyer1988)
INSTANT(boost::kreutzer1986)
INSTANT(boost::hellekalek1995)
INSTANT(boost::mt19937)
INSTANT(boost::mt11213b)
#endif
int main()
{
instantiate_all();
validate_all();
boost::mt19937 mt;
test_uniform_int(mt);
// Some compilers don't pay attention to std:3.6.1/5 and issue a
// warning here if "return 0;" is omitted.
return 0;
}

662
statistic_tests.cpp
View File

@@ -1,662 +0,0 @@
/* statistic_tests.cpp file
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
* Revision history
*/
/*
* NOTE: This is not part of the official boost submission. It exists
* only as a collection of ideas.
*/
#include <iostream>
#include <iomanip>
#include <string>
#include <functional>
#include <cmath>
#include <vector>
#include <algorithm>
#include <boost/cstdint.hpp>
#include <boost/random.hpp>
#include "statistic_tests.hpp"
#include "integrate.hpp"
namespace boost {
namespace random {
// Wikramaratna 1989 ACORN
template<class IntType, int k, IntType m, IntType val>
class additive_congruential
{
public:
typedef IntType result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const result_type min_value = 0;
static const result_type max_value = m-1;
#else
enum {
has_fixed_range = true,
min_value = 0,
max_value = m-1
};
#endif
template<class InputIterator>
explicit additive_congruential(InputIterator start) { seed(start); }
template<class InputIterator>
void seed(InputIterator start)
{
for(int i = 0; i <= k; ++i, ++start)
values[i] = *start;
}
result_type operator()()
{
for(int i = 1; i <= k; ++i) {
IntType tmp = values[i-1] + values[i];
if(tmp >= m)
tmp -= m;
values[i] = tmp;
}
return values[k];
}
result_type validation() const { return val; }
private:
IntType values[k+1];
};
template<class IntType, int r, int s, IntType m, IntType val>
class lagged_fibonacci
{
public:
typedef IntType result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const result_type min_value = 0;
static const result_type max_value = m-1;
#else
enum {
has_fixed_range = true,
min_value = 0,
max_value = m-1
};
#endif
explicit lagged_fibonacci(IntType start) { seed(start); }
template<class Generator>
explicit lagged_fibonacci(Generator & gen) { seed(gen); }
void seed(IntType start)
{
linear_congruential<uint32_t, 299375077, 0, 0, 0> init;
seed(init);
}
template<class Generator>
void seed(Generator & gen)
{
assert(r > s);
for(int i = 0; i < 607; ++i)
values[i] = gen();
current = 0;
lag = r-s;
}
result_type operator()()
{
result_type tmp = values[current] + values[lag];
if(tmp >= m)
tmp -= m;
values[current] = tmp;
++current;
if(current >= r)
current = 0;
++lag;
if(lag >= r)
lag = 0;
return tmp;
}
result_type validation() const { return val; }
private:
result_type values[r];
int current, lag;
};
} // namespace random
} // namespace boost
// distributions from Haertel's dissertation
// (additional parameterizations of the basic templates)
namespace Haertel {
typedef boost::random::linear_congruential<boost::uint64_t, 45965, 453816691,
(1ull<<31), 0> LCG_Af2;
typedef boost::random::linear_congruential<boost::uint64_t, 211936855, 0,
(1ull<<29)-3, 0> LCG_Die1;
typedef boost::random::linear_congruential<boost::uint32_t, 2824527309u, 0,
0, 0> LCG_Fis;
typedef boost::random::linear_congruential<boost::uint64_t, 950706376u, 0,
(1ull<<31)-1, 0> LCG_FM;
typedef boost::random::linear_congruential<boost::int32_t, 51081, 0,
2147483647, 0> LCG_Hae;
typedef boost::random::linear_congruential<boost::uint32_t, 69069, 1,
0, 0> LCG_VAX;
typedef boost::random::inversive_congruential<boost::int64_t, 240318, 197,
1000081, 0> NLG_Inv1;
typedef boost::random::inversive_congruential<boost::int64_t, 15707262,
13262967, (1<<24)-17, 0> NLG_Inv2;
typedef boost::random::inversive_congruential<boost::int32_t, 1, 1,
2147483647, 0> NLG_Inv4;
typedef boost::random::inversive_congruential<boost::int32_t, 1, 2,
1<<30, 0> NLG_Inv5;
typedef boost::random::additive_congruential<boost::int32_t, 6,
(1<<30)-35, 0> MRG_Acorn7;
typedef boost::random::lagged_fibonacci<boost::uint32_t, 607, 273,
0, 0> MRG_Fib2;
template<class Gen, class T>
inline void check_validation(Gen & gen, T value, const std::string & name)
{
for(int i = 0; i < 100000-1; ++i)
gen();
if(value != gen())
std::cout << name << ": validation failed" << std::endl;
}
// we have validation after 100000 steps with Haertel's generators
template<class Gen, class T>
void validate(T value, const std::string & name)
{
Gen gen(1234567);
check_validation(gen, value, name);
}
void validate_all()
{
validate<LCG_Af2>(183269031u, "LCG_Af2");
validate<LCG_Die1>(522319944u, "LCG_Die1");
validate<LCG_Fis>(-2065162233u, "LCG_Fis");
validate<LCG_FM>(581815473u, "LCG_FM");
validate<LCG_Hae>(28931709, "LCG_Hae");
validate<LCG_VAX>(1508154087u, "LCG_VAX");
validate<NLG_Inv2>(6666884, "NLG_Inv2");
validate<NLG_Inv4>(1521640076, "NLG_Inv4");
validate<NLG_Inv5>(641840839, "NLG_Inv5");
static const int acorn7_init[]
= { 1234567, 7654321, 246810, 108642, 13579, 97531, 555555 };
MRG_Acorn7 acorn7(acorn7_init);
check_validation(acorn7, 874294697, "MRG_Acorn7");
validate<MRG_Fib2>(1234567u, "MRG_Fib2");
}
} // namespace Haertel
double normal_density(double x)
{
const double pi = 3.14159265358979323846;
return 1/std::sqrt(2*pi) * std::exp(-x*x/2);
}
#ifdef _CXXRTCF_H__
namespace std {
using _CS_swamp::lgamma;
}
#endif
class chi_square_density : public std::unary_function<double, double>
{
public:
chi_square_density(int freedom)
: _exponent( static_cast<double>(freedom)/2-1 ),
_factor(1/(std::pow(2, _exponent+1) * std::exp(std::lgamma(_exponent+1))))
{ }
double operator()(double x)
{
return _factor*std::pow(x, _exponent)*std::exp(-x/2);
}
private:
double _exponent, _factor;
};
// computes F(x) or F(y) - F(x)
class chi_square_probability : public distribution_function<double>
{
public:
chi_square_probability(int freedom) : dens(freedom) {}
double operator()(double x) { return operator()(0, x); }
double operator()(double x, double y)
{ return trapezoid(dens, x, y, 1000); }
private:
chi_square_density dens;
};
class uniform_distribution : public distribution_function<double>
{
public:
uniform_distribution(double from, double to) : from(from), to(to)
{ assert(from < to); }
double operator()(double x)
{
if(x < from)
return 0;
else if(x > to)
return 1;
else
return (x-from)/(to-from);
}
double operator()(double x, double delta)
{ return operator()(x+delta) - operator()(x); }
private:
double from, to;
};
class test_environment;
class test_base
{
protected:
explicit test_base(test_environment & env) : environment(env) { }
void check(double val) const;
private:
test_environment & environment;
};
class equidistribution_test : test_base
{
public:
equidistribution_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "equidistribution: " << std::flush;
equidistribution_experiment equi(classes);
uniform_smallint<RNG> uint_linear(rng, 0, classes-1);
check(run_experiment(test_distrib_chi_square,
experiment_generator(equi, uint_linear, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(equi, uint_linear, n1), 2*n2));
std::cout << " 2D: " << std::flush;
equidistribution_2d_experiment equi_2d(classes);
unsigned int root = static_cast<unsigned int>(std::sqrt(classes));
assert(root * root == classes);
uniform_smallint<RNG> uint_square(rng, 0, root-1);
check(run_experiment(test_distrib_chi_square,
experiment_generator(equi_2d, uint_square, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(equi_2d, uint_square, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class ks_equidistribution_test : test_base
{
public:
ks_equidistribution_test(test_environment & env, unsigned int classes)
: test_base(env),
test_distrib_chi_square(kolmogorov_smirnov_probability(5000),
classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "KS: " << std::flush;
generator_reference_t<RNG> gen_ref(rng);
kolmogorov_experiment ks(n1);
uniform_distribution ud(rng.min_value, rng.max_value);
check(run_experiment(test_distrib_chi_square,
ks_experiment_generator(ks, gen_ref, ud), n2));
check(run_experiment(test_distrib_chi_square,
ks_experiment_generator(ks, gen_ref, ud), 2*n2));
}
private:
distribution_experiment test_distrib_chi_square;
};
class runs_test : test_base
{
public:
runs_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "runs: up: " << std::flush;
runs_experiment<true> r_up(classes);
generator_reference_t<RNG> gen_ref(rng);
check(run_experiment(test_distrib_chi_square,
experiment_generator(r_up, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(r_up, gen_ref, n1), 2*n2));
std::cout << " down: " << std::flush;
runs_experiment<false> r_down(classes);
check(run_experiment(test_distrib_chi_square,
experiment_generator(r_down, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(r_down, gen_ref, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class gap_test : test_base
{
public:
gap_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "gaps: " << std::flush;
gap_experiment gap(classes, 0.2, 0.8);
generator_reference_t<RNG> gen_ref(rng);
check(run_experiment(test_distrib_chi_square,
experiment_generator(gap, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(gap, gen_ref, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class poker_test : test_base
{
public:
poker_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "poker: " << std::flush;
poker_experiment poker(8, classes);
uniform_smallint<RNG> usmall(rng, 0, 7);
check(run_experiment(test_distrib_chi_square,
experiment_generator(poker, usmall, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(poker, usmall, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class coupon_collector_test : test_base
{
public:
coupon_collector_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "coupon collector: " << std::flush;
coupon_collector_experiment coupon(5, classes);
uniform_smallint<RNG> usmall(rng, 0, 4);
check(run_experiment(test_distrib_chi_square,
experiment_generator(coupon, usmall, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(coupon, usmall, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class permutation_test : test_base
{
public:
permutation_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(fac<int>(classes)-1),
high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "permutation: " << std::flush;
permutation_experiment perm(classes);
generator_reference_t<RNG> gen_ref(rng);
check(run_experiment(test_distrib_chi_square,
experiment_generator(perm, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(perm, gen_ref, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class maximum_test : test_base
{
public:
maximum_test(test_environment & env, unsigned int high_classes)
: test_base(env),
test_distrib_chi_square(kolmogorov_smirnov_probability(1000),
high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "maximum-of-t: " << std::flush;
maximum_experiment<RNG> mx(rng, n1, 5);
check(run_experiment(test_distrib_chi_square, mx, n2));
check(run_experiment(test_distrib_chi_square, mx, 2*n2));
std::cout << std::endl;
}
private:
distribution_experiment test_distrib_chi_square;
};
class birthday_test : test_base
{
public:
birthday_test(test_environment & env)
: test_base(env)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "birthday spacing: " << std::flush;
uniform_int<RNG> uni(rng, 0, (1<<25));
birthday_spacing_experiment bsp(4, 512, (1<<25));
std::cout << run_experiment(bsp, uni, n1);
#if 0
check(run_experiment(test_distrib_chi_square,
experiment_generator(perm, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(perm, gen_ref, n1), 2*n2));
#endif
std::cout << std::endl;
}
};
class test_environment
{
public:
static const int classes = 20;
explicit test_environment(double confid)
: confidence(confid),
confidence_chi_square_quantil(quantil(chi_square_density(classes-1), 0, confidence, 1e-4)),
test_distrib_chi_square6(chi_square_probability(7-1), classes),
ksequi_test(*this, classes),
equi_test(*this, 100, classes),
rns_test(*this, 7, classes),
gp_test(*this, 7, classes),
pk_test(*this, 5, classes),
cpn_test(*this, 15, classes),
perm_test(*this, 5, classes),
max_test(*this, classes),
bday_test(*this)
{
std::cout << "Confidence level: " << confid
<< "; 1-alpha = " << (1-confid)
<< "; chi_square(" << (classes-1)
<< ", " << confidence_chi_square_quantil
<< ") = "
<< chi_square_probability(classes-1)(0, confidence_chi_square_quantil)
<< std::endl;
}
bool check_confidence(double val, double chi_square_conf) const
{
std::cout << val;
bool result = (val <= chi_square_conf);
if(!result) {
std::cout << "* [";
double prob = (val > 10*chi_square_conf ? 1 :
chi_square_probability(classes-1)(0, val));
std::cout << (1-prob) << "]";
}
std::cout << " " << std::flush;
return result;
}
bool check(double chi_square_value) const
{
return check_confidence(chi_square_value, confidence_chi_square_quantil);
}
template<class RNG>
void run_test(const std::string & name)
{
using namespace boost;
std::cout << "Running tests on " << name << std::endl;
RNG rng(1234567);
typedef boost::uniform_01<RNG> UGen;
#if 1
ksequi_test.run(rng, 5000, 250);
equi_test.run(rng, 5000, 250);
rns_test.run(rng, 100000, 250);
gp_test.run(rng, 10000, 250);
pk_test.run(rng, 5000, 250);
cpn_test.run(rng, 500, 250);
perm_test.run(rng, 1200, 250);
max_test.run(rng, 1000, 250);
#endif
bday_test.run(rng, 1000, 150);
std::cout << std::endl;
}
private:
double confidence;
double confidence_chi_square_quantil;
distribution_experiment test_distrib_chi_square6;
ks_equidistribution_test ksequi_test;
equidistribution_test equi_test;
runs_test rns_test;
gap_test gp_test;
poker_test pk_test;
coupon_collector_test cpn_test;
permutation_test perm_test;
maximum_test max_test;
birthday_test bday_test;
};
void test_base::check(double val) const
{
environment.check(val);
}
void print_ks_table()
{
std::cout.setf(std::ios::fixed);
std::cout.precision(5);
static const double all_p[] = { 0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99 };
for(int n = 0; n <= 10000; (n < 55 ? ++n : n *= 10)) {
std::cout << std::setw(4) << n << " ";
for(unsigned int i = 0; i < sizeof(all_p)/sizeof(all_p[0]); ++i) {
std::cout << std::setw(8)
<< (n == 0 ? all_p[i] :
invert_monotone_inc(kolmogorov_smirnov_probability(n), all_p[i], 0, 10))
<< " ";
}
std::cout << std::endl;
}
}
int main()
{
// Haertel::validate_all();
test_environment env(0.99);
env.run_test<boost::minstd_rand>("minstd_rand");
env.run_test<boost::mt19937>("mt19937");
env.run_test<Haertel::LCG_Af2>("LCG_Af2");
env.run_test<Haertel::LCG_Die1>("LCG_Die1");
env.run_test<Haertel::LCG_Fis>("LCG_Fis");
env.run_test<Haertel::LCG_FM>("LCG_FM");
env.run_test<Haertel::LCG_Hae>("LCG_Hae");
env.run_test<Haertel::LCG_VAX>("LCG_VAX");
env.run_test<Haertel::NLG_Inv1>("NLG_Inv1");
env.run_test<Haertel::NLG_Inv2>("NLG_Inv2");
env.run_test<Haertel::NLG_Inv4>("NLG_Inv4");
env.run_test<Haertel::NLG_Inv5>("NLG_Inv5");
}

628
statistic_tests.hpp
View File

@@ -1,628 +0,0 @@
/* statistic_tests.hpp header file
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without free provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
*/
#ifndef STATISTIC_TESTS_HPP
#define STATISTIC_TESTS_HPP
#include <stdexcept>
#include <iterator>
#include <vector>
#include <limits>
#include <algorithm>
#include <boost/random.hpp>
template<class T>
inline T fac(int k)
{
T result = 1;
for(T i = 2; i <= k; ++i)
result *= i;
return result;
}
template<class T>
T binomial(int n, int k)
{
if(k < n/2)
k = n-k;
T result = 1;
for(int i = k+1; i<= n; ++i)
result *= i;
return result / fac<T>(n-k);
}
template<class T>
T stirling2(int n, int m)
{
T sum = 0;
for(int k = 0; k <= m; ++k)
sum += binomial<T>(m, k) * std::pow(k, n) * ( (m-k)%2 == 0 ? 1 : -1);
return sum / fac<T>(m);
}
/*
* Experiments which create an empirical distribution in classes,
* suitable for the chi-square test.
*/
// std::floor(gen() * classes)
class experiment_base
{
public:
experiment_base(int cls) : _classes(cls) { }
unsigned int classes() const { return _classes; }
protected:
unsigned int _classes;
};
class equidistribution_experiment : public experiment_base
{
public:
explicit equidistribution_experiment(unsigned int classes)
: experiment_base(classes) { }
template<class NumberGenerator, class Counter>
void run(NumberGenerator f, Counter & count, int n) const
{
assert(f.min() == 0 &&
static_cast<unsigned int>(f.max()) == classes()-1);
for(int i = 0; i < n; ++i)
count(f());
}
double probability(int i) const { return 1.0/classes(); }
};
// two-dimensional equidistribution experiment
class equidistribution_2d_experiment : public equidistribution_experiment
{
public:
explicit equidistribution_2d_experiment(unsigned int classes)
: equidistribution_experiment(classes) { }
template<class NumberGenerator, class Counter>
void run(NumberGenerator f, Counter & count, int n) const
{
unsigned int range = f.max()+1;
assert(f.min() == 0 && range*range == classes());
for(int i = 0; i < n; ++i) {
int y1 = f();
int y2 = f();
count(y1 + range * y2);
}
}
};
// distribution experiment: assume a probability density and
// count events so that an equidistribution results.
class distribution_experiment : public equidistribution_experiment
{
public:
template<class UnaryFunction>
distribution_experiment(UnaryFunction probability , unsigned int classes)
: equidistribution_experiment(classes), limit(classes)
{
for(unsigned int i = 0; i < classes-1; ++i)
limit[i] = invert_monotone_inc(probability, (i+1)*0.05, 0, 1000);
limit[classes-1] = std::numeric_limits<double>::infinity();
#if 0
std::cout << __PRETTY_FUNCTION__ << ": ";
for(unsigned int i = 0; i < classes; ++i)
std::cout << limit[i] << " ";
std::cout << std::endl;
#endif
}
template<class NumberGenerator, class Counter>
void run(NumberGenerator f, Counter & count, int n) const
{
for(int i = 0; i < n; ++i) {
limits_type::const_iterator it =
std::lower_bound(limit.begin(), limit.end(), f());
count(it-limit.begin());
}
}
private:
typedef std::vector<double> limits_type;
limits_type limit;
};
// runs-up/runs-down experiment
template<bool up>
class runs_experiment : public experiment_base
{
public:
explicit runs_experiment(unsigned int classes) : experiment_base(classes) { }
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
result_type init = (up ? f.min() : f.max());
result_type previous = init;
unsigned int length = 0;
for(int i = 0; i < n; ++i) {
result_type val = f();
if(up ? previous <= val : previous >= val) {
previous = val;
++length;
} else {
count(std::min(length, classes())-1);
length = 0;
previous = init;
// don't use this value, so that runs are independent
}
}
}
double probability(unsigned int r) const
{
if(r == classes()-1)
return 1.0/fac<double>(classes());
else
return static_cast<double>(r+1)/fac<double>(r+2);
}
};
// gap length experiment
class gap_experiment : public experiment_base
{
public:
gap_experiment(unsigned int classes, double alpha, double beta)
: experiment_base(classes), alpha(alpha), beta(beta) { }
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
double range = f.max() - f.min() + 1.0;
result_type low = static_cast<result_type>(alpha * range);
result_type high = static_cast<result_type>(beta * range);
unsigned int length = 0;
for(int i = 0; i < n; ) {
result_type value = f() - f.min();
if(value < low || value > high)
++length;
else {
count(std::min(length, classes()-1));
length = 0;
++i;
}
}
}
double probability(unsigned int r) const
{
double p = beta-alpha;
if(r == classes()-1)
return std::pow(1-p, static_cast<double>(r));
else
return p * std::pow(1-p, static_cast<double>(r));
}
private:
double alpha, beta;
};
// poker experiment
class poker_experiment : public experiment_base
{
public:
poker_experiment(unsigned int d, unsigned int k)
: experiment_base(k), range(d)
{
assert(range > 1);
}
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
assert(std::numeric_limits<result_type>::is_integer);
assert(f.min() == 0);
assert(f.max() == static_cast<result_type>(range-1));
std::vector<result_type> v(classes());
for(int i = 0; i < n; ++i) {
for(unsigned int j = 0; j < classes(); ++j)
v[j] = f();
std::sort(v.begin(), v.end());
result_type prev = v[0];
int r = 1; // count different values in v
for(unsigned int i = 1; i < classes(); ++i) {
if(prev != v[i]) {
prev = v[i];
++r;
}
}
count(r-1);
}
}
double probability(unsigned int r) const
{
++r; // transform to 1 <= r <= 5
double result = range;
for(unsigned int i = 1; i < r; ++i)
result *= range-i;
return result / std::pow(range, static_cast<double>(classes())) *
stirling2<double>(classes(), r);
}
private:
unsigned int range;
};
// coupon collector experiment
class coupon_collector_experiment : public experiment_base
{
public:
coupon_collector_experiment(unsigned int d, unsigned int cls)
: experiment_base(cls), d(d)
{
assert(d > 1);
}
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
assert(std::numeric_limits<result_type>::is_integer);
assert(f.min() == 0);
assert(f.max() == static_cast<result_type>(d-1));
std::vector<bool> occurs(d);
for(int i = 0; i < n; ++i) {
occurs.assign(d, false);
unsigned int r = 0; // length of current sequence
int q = 0; // number of non-duplicates in current set
for(;;) {
result_type val = f();
++r;
if(!occurs[val]) { // new set element
occurs[val] = true;
++q;
if(q == d)
break; // one complete set
}
}
count(std::min(r-d, classes()-1));
}
}
double probability(unsigned int r) const
{
if(r == classes()-1)
return 1-fac<double>(d)/std::pow(d, static_cast<double>(d+classes()-2))*
stirling2<double>(d+classes()-2, d);
else
return fac<double>(d)/std::pow(d, static_cast<double>(d+r)) *
stirling2<double>(d+r-1, d-1);
}
private:
int d;
};
// permutation test
class permutation_experiment : public equidistribution_experiment
{
public:
permutation_experiment(unsigned int t)
: equidistribution_experiment(fac<int>(t)), t(t)
{
assert(t > 1);
}
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
std::vector<result_type> v(t);
for(int i = 0; i < n; ++i) {
std::generate_n(v.begin(), t, f);
int x = 0;
for(int r = t-1; r > 0; r--) {
std::vector<result_type>::iterator it =
std::max_element(v.begin(), v.begin()+r+1);
x = (r+1)*x + (it-v.begin());
std::iter_swap(it, v.begin()+r);
}
count(x);
}
}
private:
int t;
};
// birthday spacing experiment test
class birthday_spacing_experiment : public experiment_base
{
public:
birthday_spacing_experiment(unsigned int d, int n, int m)
: experiment_base(d), n(n), m(m)
{
}
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n_total) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
assert(std::numeric_limits<result_type>::is_integer);
assert(f.min() == 0);
assert(f.max() == static_cast<result_type>(m-1));
for(int j = 0; j < n_total; j++) {
std::vector<result_type> v(n);
std::generate_n(v.begin(), n, f);
std::sort(v.begin(), v.end());
std::vector<result_type> spacing(n);
for(int i = 0; i < n-1; i++)
spacing[i] = v[i+1]-v[i];
spacing[n-1] = v[0] + m - v[n-1];
std::sort(spacing.begin(), spacing.end());
unsigned int k = 0;
for(int i = 0; i < n-1; ++i) {
if(spacing[i] == spacing[i+1])
++k;
}
count(std::min(k, classes()-1));
}
}
double probability(unsigned int r) const
{
assert(classes() == 4);
assert(m == (1<<25));
assert(n == 512);
static const double prob[] = { 0.368801577, 0.369035243, 0.183471182,
0.078691997 };
return prob[r];
}
private:
int n, m;
};
/*
* Misc. helper functions.
*/
template<class Float>
struct distribution_function
{
typedef Float result_type;
typedef Float argument_type;
typedef Float first_argument_type;
typedef Float second_argument_type;
};
// computes P(K_n <= t) or P(t1 <= K_n <= t2). See Knuth, 3.3.1
class kolmogorov_smirnov_probability : public distribution_function<double>
{
public:
kolmogorov_smirnov_probability(int n)
: approx(n > 50), n(n), sqrt_n(std::sqrt(n))
{
if(!approx)
n_n = std::pow(static_cast<double>(n), n);
}
double operator()(double t) const
{
if(approx) {
return 1-std::exp(-2*t*t)*(1-2.0/3.0*t/sqrt_n);
} else {
t *= sqrt_n;
double sum = 0;
for(int k = static_cast<int>(std::ceil(t)); k <= n; k++)
sum += binomial<double>(n, k) * std::pow(k-t, k) *
std::pow(t+n-k, n-k-1);
return 1 - t/n_n * sum;
}
}
double operator()(double t1, double t2) const
{ return operator()(t2) - operator()(t1); }
private:
bool approx;
int n;
double sqrt_n;
double n_n;
};
/*
* Experiments for generators with continuous distribution functions
*/
class kolmogorov_experiment
{
public:
kolmogorov_experiment(int n) : n(n), ksp(n) { }
template<class NumberGenerator, class Distribution>
double run(NumberGenerator gen, Distribution distrib) const
{
const int m = n;
typedef std::vector<double> saved_temp;
saved_temp a(m,1.0), b(m,0);
std::vector<int> c(m,0);
for(int i = 0; i < n; ++i) {
double val = gen();
double y = distrib(val);
int k = static_cast<int>(std::floor(m*y));
if(k >= m)
--k; // should not happen
a[k] = std::min(a[k], y);
b[k] = std::max(b[k], y);
++c[k];
}
double kplus = 0, kminus = 0;
int j = 0;
for(int k = 0; k < m; ++k) {
if(c[k] > 0) {
kminus = std::max(kminus, a[k]-j/static_cast<double>(n));
j += c[k];
kplus = std::max(kplus, j/static_cast<double>(n) - b[k]);
}
}
kplus *= std::sqrt(n);
kminus *= std::sqrt(n);
// std::cout << "k+ " << kplus << " k- " << kminus << std::endl;
return kplus;
}
double probability(double x) const
{
return ksp(x);
}
private:
int n;
kolmogorov_smirnov_probability ksp;
};
// maximum-of-t test (KS-based)
template<class UniformRandomNumberGenerator>
class maximum_experiment
{
public:
typedef UniformRandomNumberGenerator base_type;
maximum_experiment(base_type & f, int n, int t) : f(f), ke(n), t(t)
{ }
double operator()() const
{
double res = ke.run(generator<base_type>(f, t),
std::bind2nd(std::ptr_fun(static_cast<double (*)(double, double)>(&std::pow)), t));
return res;
}
private:
struct generator {
generator(base_type & f, int t) : f(f), t(t) { }
double operator()()
{
double mx = f();
for(int i = 1; i < t; ++i)
mx = std::max(mx, f());
return mx;
}
private:
boost::uniform_01<base_type> f;
int t;
};
base_type & f;
kolmogorov_experiment ke;
int t;
};
// compute a chi-square value for the distribution approximation error
template<class ForwardIterator, class UnaryFunction>
typename UnaryFunction::result_type
chi_square_value(ForwardIterator first, ForwardIterator last,
UnaryFunction probability)
{
typedef std::iterator_traits<ForwardIterator> iter_traits;
typedef typename iter_traits::value_type counter_type;
typedef typename UnaryFunction::result_type result_type;
unsigned int classes = std::distance(first, last);
result_type sum = 0;
counter_type n = 0;
for(unsigned int i = 0; i < classes; ++first, ++i) {
counter_type count = *first;
n += count;
sum += (count/probability(i)) * count; // avoid overflow
}
#if 0
for(unsigned int i = 0; i < classes; ++i) {
// std::cout << (n*probability(i)) << " ";
if(n * probability(i) < 5)
std::cerr << "Not enough test runs for slot " << i
<< " p=" << probability(i) << ", n=" << n
<< std::endl;
}
#endif
// std::cout << std::endl;
// throw std::invalid_argument("not enough test runs");
return sum/n - n;
}
template<class RandomAccessContainer>
class generic_counter
{
public:
explicit generic_counter(unsigned int classes) : container(classes, 0) { }
void operator()(int i)
{
assert(i >= 0);
assert(static_cast<unsigned int>(i) < container.size());
++container[i];
}
typename RandomAccessContainer::const_iterator begin() const
{ return container.begin(); }
typename RandomAccessContainer::const_iterator end() const
{ return container.end(); }
private:
RandomAccessContainer container;
};
// chi_square test
template<class Experiment, class Generator>
double run_experiment(const Experiment & experiment, Generator gen, int n)
{
generic_counter<std::vector<int> > v(experiment.classes());
experiment.run(gen, v, n);
return chi_square_value(v.begin(), v.end(),
std::bind1st(std::mem_fun_ref(&Experiment::probability),
experiment));
}
// number generator with experiment results (for nesting)
template<class Experiment, class Generator>
class experiment_generator_t
{
public:
experiment_generator_t(const Experiment & exper, Generator & gen, int n)
: experiment(exper), generator(gen), n(n) { }
double operator()() { return run_experiment(experiment, generator, n); }
private:
const Experiment & experiment;
Generator & generator;
int n;
};
template<class Experiment, class Generator>
experiment_generator_t<Experiment, Generator>
experiment_generator(const Experiment & e, Generator & gen, int n)
{
return experiment_generator_t<Experiment, Generator>(e, gen, n);
}
template<class Experiment, class Generator, class Distribution>
class ks_experiment_generator_t
{
public:
ks_experiment_generator_t(const Experiment & exper, Generator & gen,
const Distribution & distrib)
: experiment(exper), generator(gen), distribution(distrib) { }
double operator()() { return experiment.run(generator, distribution); }
private:
const Experiment & experiment;
Generator & generator;
Distribution distribution;
};
template<class Experiment, class Generator, class Distribution>
ks_experiment_generator_t<Experiment, Generator, Distribution>
ks_experiment_generator(const Experiment & e, Generator & gen,
const Distribution & distrib)
{
return ks_experiment_generator_t<Experiment, Generator, Distribution>
(e, gen, distrib);
}
#endif /* STATISTIC_TESTS_HPP */