Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: downloads/boost_1_34_1/libs/random/statistic_tests.hpp @ 29

Last change on this file since 29 was 29, checked in by landauf, 16 years ago

updated boost from 1_33_1 to 1_34_1

File size: 17.3 KB
Line 
1/* statistic_tests.hpp header file
2 *
3 * Copyright Jens Maurer 2000
4 * Distributed under the Boost Software License, Version 1.0. (See
5 * accompanying file LICENSE_1_0.txt or copy at
6 * http://www.boost.org/LICENSE_1_0.txt)
7 *
8 * $Id: statistic_tests.hpp,v 1.13 2004/07/27 03:43:34 dgregor Exp $
9 *
10 */
11
12#ifndef STATISTIC_TESTS_HPP
13#define STATISTIC_TESTS_HPP
14
15#include <stdexcept>
16#include <iterator>
17#include <vector>
18#include <boost/limits.hpp>
19#include <algorithm>
20#include <cmath>
21
22#include <boost/random.hpp>
23#include <boost/config.hpp>
24
25
26#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300
27namespace std
28{
29  inline double pow(double a, double b) { return ::pow(a,b); }
30  inline double ceil(double x) { return ::ceil(x); }
31} // namespace std
32#endif
33
34
35template<class T>
36inline T fac(int k)
37{
38  T result = 1;
39  for(T i = 2; i <= k; ++i)
40    result *= i;
41  return result;
42}
43
44template<class T>
45T binomial(int n, int k)
46{
47  if(k < n/2)
48    k = n-k;
49  T result = 1;
50  for(int i = k+1; i<= n; ++i)
51    result *= i;
52  return result / fac<T>(n-k);
53}
54
55template<class T>
56T stirling2(int n, int m)
57{
58  T sum = 0;
59  for(int k = 0; k <= m; ++k)
60    sum += binomial<T>(m, k) * std::pow(double(k), n) *
61      ( (m-k)%2 == 0 ? 1 : -1);
62  return sum / fac<T>(m);
63}
64
65/*
66 * Experiments which create an empirical distribution in classes,
67 * suitable for the chi-square test.
68 */
69// std::floor(gen() * classes)
70
71class experiment_base
72{
73public:
74  experiment_base(int cls) : _classes(cls) { }
75  unsigned int classes() const { return _classes; }
76protected:
77  unsigned int _classes;
78};
79
80class equidistribution_experiment : public experiment_base
81{
82public:
83  explicit equidistribution_experiment(unsigned int classes) 
84    : experiment_base(classes) { }
85 
86  template<class NumberGenerator, class Counter>
87  void run(NumberGenerator f, Counter & count, int n) const
88  {
89    assert((f.min)() == 0 &&
90           static_cast<unsigned int>((f.max)()) == classes()-1);
91    for(int i = 0; i < n; ++i)
92      count(f());
93  }
94  double probability(int i) const { return 1.0/classes(); }
95};
96
97// two-dimensional equidistribution experiment
98class equidistribution_2d_experiment : public equidistribution_experiment
99{
100public:
101  explicit equidistribution_2d_experiment(unsigned int classes) 
102    : equidistribution_experiment(classes) { }
103
104  template<class NumberGenerator, class Counter>
105  void run(NumberGenerator f, Counter & count, int n) const
106  {
107    unsigned int range = (f.max)()+1;
108    assert((f.min)() == 0 && range*range == classes());
109    for(int i = 0; i < n; ++i) {
110      int y1 = f();
111      int y2 = f();
112      count(y1 + range * y2);
113    }
114  }
115};
116
117// distribution experiment: assume a probability density and
118// count events so that an equidistribution results.
119class distribution_experiment : public equidistribution_experiment
120{
121public:
122  template<class UnaryFunction>
123  distribution_experiment(UnaryFunction probability , unsigned int classes)
124    : equidistribution_experiment(classes), limit(classes)
125  {
126    for(unsigned int i = 0; i < classes-1; ++i)
127      limit[i] = invert_monotone_inc(probability, (i+1)*0.05, 0, 1000);
128    limit[classes-1] = std::numeric_limits<double>::infinity();
129    if(limit[classes-1] < (std::numeric_limits<double>::max)())
130      limit[classes-1] = (std::numeric_limits<double>::max)();
131#if 0
132    std::cout << __PRETTY_FUNCTION__ << ": ";
133    for(unsigned int i = 0; i < classes; ++i)
134      std::cout << limit[i] << " ";
135    std::cout << std::endl;
136#endif
137  }
138
139  template<class NumberGenerator, class Counter>
140  void run(NumberGenerator f, Counter & count, int n) const
141  {
142    for(int i = 0; i < n; ++i) {
143      limits_type::const_iterator it =
144        std::lower_bound(limit.begin(), limit.end(), f());
145      count(it-limit.begin());
146    }
147  }
148private:
149  typedef std::vector<double> limits_type;
150  limits_type limit;
151};
152
153// runs-up/runs-down experiment
154template<bool up>
155class runs_experiment : public experiment_base
156{
157public:
158  explicit runs_experiment(unsigned int classes) : experiment_base(classes) { }
159 
160  template<class UniformRandomNumberGenerator, class Counter>
161  void run(UniformRandomNumberGenerator f, Counter & count, int n) const
162  {
163    typedef typename UniformRandomNumberGenerator::result_type result_type;
164    result_type init = (up ? (f.min)() : (f.max)());
165    result_type previous = init;
166    unsigned int length = 0;
167    for(int i = 0; i < n; ++i) {
168      result_type val = f();
169      if(up ? previous <= val : previous >= val) {
170        previous = val;
171        ++length;
172      } else {
173        count((std::min)(length, classes())-1);
174        length = 0;
175        previous = init;
176        // don't use this value, so that runs are independent
177      }
178    }
179  }
180  double probability(unsigned int r) const
181  {
182    if(r == classes()-1)
183      return 1.0/fac<double>(classes());
184    else
185      return static_cast<double>(r+1)/fac<double>(r+2);
186  }
187};
188
189// gap length experiment
190class gap_experiment : public experiment_base
191{
192public:
193  gap_experiment(unsigned int classes, double alpha, double beta)
194    : experiment_base(classes), alpha(alpha), beta(beta) { }
195 
196  template<class UniformRandomNumberGenerator, class Counter>
197  void run(UniformRandomNumberGenerator f, Counter & count, int n) const
198  {
199    typedef typename UniformRandomNumberGenerator::result_type result_type;
200    double range = (f.max)() - (f.min)() + 1.0;
201    result_type low = static_cast<result_type>(alpha * range);
202    result_type high = static_cast<result_type>(beta * range);
203    unsigned int length = 0;
204    for(int i = 0; i < n; ) {
205      result_type value = f() - (f.min)();
206      if(value < low || value > high)
207        ++length;
208      else {
209        count((std::min)(length, classes()-1));
210        length = 0;
211        ++i;
212      }
213    }
214  }
215  double probability(unsigned int r) const
216  {
217    double p = beta-alpha;
218    if(r == classes()-1)
219      return std::pow(1-p, static_cast<double>(r));
220    else
221      return p * std::pow(1-p, static_cast<double>(r));
222  }
223private:
224  double alpha, beta;
225};
226
227// poker experiment
228class poker_experiment : public experiment_base
229{
230public:
231  poker_experiment(unsigned int d, unsigned int k)
232    : experiment_base(k), range(d)
233  {
234    assert(range > 1);
235  }
236
237  template<class UniformRandomNumberGenerator, class Counter>
238  void run(UniformRandomNumberGenerator f, Counter & count, int n) const
239  {
240    typedef typename UniformRandomNumberGenerator::result_type result_type;
241    assert(std::numeric_limits<result_type>::is_integer);
242    assert((f.min)() == 0);
243    assert((f.max)() == static_cast<result_type>(range-1));
244    std::vector<result_type> v(classes());
245    for(int i = 0; i < n; ++i) {
246      for(unsigned int j = 0; j < classes(); ++j)
247        v[j] = f();
248      std::sort(v.begin(), v.end());
249      result_type prev = v[0];
250      int r = 1;     // count different values in v
251      for(unsigned int i = 1; i < classes(); ++i) {
252        if(prev != v[i]) {
253          prev = v[i];
254          ++r;
255        }
256      }
257      count(r-1);
258    }
259  }
260
261  double probability(unsigned int r) const
262  {
263    ++r;       // transform to 1 <= r <= 5
264    double result = range;
265    for(unsigned int i = 1; i < r; ++i)
266      result *= range-i;
267    return result / std::pow(range, static_cast<double>(classes())) *
268      stirling2<double>(classes(), r);
269  }
270private:
271  unsigned int range;
272};
273
274// coupon collector experiment
275class coupon_collector_experiment : public experiment_base
276{
277public:
278  coupon_collector_experiment(unsigned int d, unsigned int cls)
279    : experiment_base(cls), d(d)
280  {
281    assert(d > 1);
282  }
283
284  template<class UniformRandomNumberGenerator, class Counter>
285  void run(UniformRandomNumberGenerator f, Counter & count, int n) const
286  {
287    typedef typename UniformRandomNumberGenerator::result_type result_type;
288    assert(std::numeric_limits<result_type>::is_integer);
289    assert((f.min)() == 0);
290    assert((f.max)() == static_cast<result_type>(d-1));
291    std::vector<bool> occurs(d);
292    for(int i = 0; i < n; ++i) {
293      occurs.assign(d, false);
294      unsigned int r = 0;            // length of current sequence
295      int q = 0;                     // number of non-duplicates in current set
296      for(;;) {
297        result_type val = f();
298        ++r;
299        if(!occurs[val]) {       // new set element
300          occurs[val] = true;
301          ++q;
302          if(q == d)
303            break;     // one complete set
304        }
305      }
306      count((std::min)(r-d, classes()-1));
307    }
308  }
309  double probability(unsigned int r) const
310  {
311    if(r == classes()-1)
312      return 1-fac<double>(d)/std::pow(d, static_cast<double>(d+classes()-2))* 
313        stirling2<double>(d+classes()-2, d);
314    else
315      return fac<double>(d)/std::pow(d, static_cast<double>(d+r)) * 
316        stirling2<double>(d+r-1, d-1);
317  }
318private:
319  int d;
320};
321
322// permutation test
323class permutation_experiment : public equidistribution_experiment
324{
325public:
326  permutation_experiment(unsigned int t)
327    : equidistribution_experiment(fac<int>(t)), t(t)
328  {
329    assert(t > 1);
330  }
331
332  template<class UniformRandomNumberGenerator, class Counter>
333  void run(UniformRandomNumberGenerator f, Counter & count, int n) const
334  {
335    typedef typename UniformRandomNumberGenerator::result_type result_type;
336    std::vector<result_type> v(t);
337    for(int i = 0; i < n; ++i) {
338      std::generate_n(v.begin(), t, f);
339      int x = 0;
340      for(int r = t-1; r > 0; r--) {
341        typename std::vector<result_type>::iterator it = 
342          std::max_element(v.begin(), v.begin()+r+1);
343        x = (r+1)*x + (it-v.begin());
344        std::iter_swap(it, v.begin()+r);
345      }
346      count(x);
347    }
348  }
349private:
350  int t;
351};
352
353// birthday spacing experiment test
354class birthday_spacing_experiment : public experiment_base
355{
356public:
357  birthday_spacing_experiment(unsigned int d, int n, int m)
358    : experiment_base(d), n(n), m(m)
359  {
360  }
361
362  template<class UniformRandomNumberGenerator, class Counter>
363  void run(UniformRandomNumberGenerator f, Counter & count, int n_total) const
364  {
365    typedef typename UniformRandomNumberGenerator::result_type result_type;
366    assert(std::numeric_limits<result_type>::is_integer);
367    assert((f.min)() == 0);
368    assert((f.max)() == static_cast<result_type>(m-1));
369   
370    for(int j = 0; j < n_total; j++) {
371      std::vector<result_type> v(n);
372      std::generate_n(v.begin(), n, f);
373      std::sort(v.begin(), v.end());
374      std::vector<result_type> spacing(n);
375      for(int i = 0; i < n-1; i++)
376        spacing[i] = v[i+1]-v[i];
377      spacing[n-1] = v[0] + m - v[n-1];
378      std::sort(spacing.begin(), spacing.end());
379      unsigned int k = 0;
380      for(int i = 0; i < n-1; ++i) {
381        if(spacing[i] == spacing[i+1])
382          ++k;
383      }
384      count((std::min)(k, classes()-1));
385    }
386  }
387
388  double probability(unsigned int r) const
389  {
390    assert(classes() == 4);
391    assert(m == (1<<25));
392    assert(n == 512);
393    static const double prob[] = { 0.368801577, 0.369035243, 0.183471182,
394                                   0.078691997 };
395    return prob[r];
396  }
397private:
398  int n, m;
399};
400/*
401 * Misc. helper functions.
402 */
403
404template<class Float>
405struct distribution_function
406{
407  typedef Float result_type;
408  typedef Float argument_type;
409  typedef Float first_argument_type;
410  typedef Float second_argument_type;
411};
412
413// computes P(K_n <= t) or P(t1 <= K_n <= t2).  See Knuth, 3.3.1
414class kolmogorov_smirnov_probability : public distribution_function<double>
415{
416public:
417  kolmogorov_smirnov_probability(int n) 
418    : approx(n > 50), n(n), sqrt_n(std::sqrt(double(n)))
419  {
420    if(!approx)
421      n_n = std::pow(static_cast<double>(n), n);
422  }
423 
424  double operator()(double t) const
425  {
426    if(approx) {
427      return 1-std::exp(-2*t*t)*(1-2.0/3.0*t/sqrt_n);
428    } else {
429      t *= sqrt_n;
430      double sum = 0;
431      for(int k = static_cast<int>(std::ceil(t)); k <= n; k++)
432        sum += binomial<double>(n, k) * std::pow(k-t, k) * 
433          std::pow(t+n-k, n-k-1);
434      return 1 - t/n_n * sum;
435    }
436  }
437  double operator()(double t1, double t2) const
438  { return operator()(t2) - operator()(t1); }
439
440private:
441  bool approx;
442  int n;
443  double sqrt_n;
444  double n_n;
445};
446
447/*
448 * Experiments for generators with continuous distribution functions
449 */
450class kolmogorov_experiment
451{
452public:
453  kolmogorov_experiment(int n) : n(n), ksp(n) { }
454  template<class NumberGenerator, class Distribution>
455  double run(NumberGenerator gen, Distribution distrib) const
456  {
457    const int m = n;
458    typedef std::vector<double> saved_temp;
459    saved_temp a(m,1.0), b(m,0);
460    std::vector<int> c(m,0);
461    for(int i = 0; i < n; ++i) {
462      double val = gen();
463      double y = distrib(val);
464      int k = static_cast<int>(std::floor(m*y));
465      if(k >= m)
466        --k;    // should not happen
467      a[k] = (std::min)(a[k], y);
468      b[k] = (std::max)(b[k], y);
469      ++c[k];
470    }
471    double kplus = 0, kminus = 0;
472    int j = 0;
473    for(int k = 0; k < m; ++k) {
474      if(c[k] > 0) {
475        kminus = (std::max)(kminus, a[k]-j/static_cast<double>(n));
476        j += c[k];
477        kplus = (std::max)(kplus, j/static_cast<double>(n) - b[k]);
478      }
479    }
480    kplus *= std::sqrt(double(n));
481    kminus *= std::sqrt(double(n));
482    // std::cout << "k+ " << kplus << "   k- " << kminus << std::endl;
483    return kplus;
484  }
485  double probability(double x) const
486  {
487    return ksp(x);
488  }
489private:
490  int n;
491  kolmogorov_smirnov_probability ksp;
492};
493
494// maximum-of-t test (KS-based)
495template<class UniformRandomNumberGenerator>
496class maximum_experiment
497{
498public:
499  typedef UniformRandomNumberGenerator base_type;
500  maximum_experiment(base_type & f, int n, int t) : f(f), ke(n), t(t)
501  { }
502
503  double operator()() const
504  {
505    double res = ke.run(generator(f, t), 
506                  std::bind2nd(std::ptr_fun(static_cast<double (*)(double, double)>(&std::pow)), t));
507    return res;
508  }
509
510private:
511  struct generator {
512    generator(base_type & f, int t) : f(f), t(t) { }
513    double operator()()
514    {
515      double mx = f();
516      for(int i = 1; i < t; ++i)
517        mx = (std::max)(mx, f());
518      return mx;
519    }
520  private:
521    boost::uniform_01<base_type> f;
522    int t;
523  };
524  base_type & f;
525  kolmogorov_experiment ke;
526  int t;
527};
528
529// compute a chi-square value for the distribution approximation error
530template<class ForwardIterator, class UnaryFunction>
531typename UnaryFunction::result_type
532chi_square_value(ForwardIterator first, ForwardIterator last,
533                 UnaryFunction probability)
534{
535  typedef std::iterator_traits<ForwardIterator> iter_traits;
536  typedef typename iter_traits::value_type counter_type;
537  typedef typename UnaryFunction::result_type result_type;
538  unsigned int classes = std::distance(first, last);
539  result_type sum = 0;
540  counter_type n = 0;
541  for(unsigned int i = 0; i < classes; ++first, ++i) {
542    counter_type count = *first;
543    n += count;
544    sum += (count/probability(i)) * count;  // avoid overflow
545  }
546#if 0
547  for(unsigned int i = 0; i < classes; ++i) {
548    // std::cout << (n*probability(i)) << " ";
549    if(n * probability(i) < 5)
550      std::cerr << "Not enough test runs for slot " << i
551                << " p=" << probability(i) << ", n=" << n
552                << std::endl;
553  }
554#endif
555  // std::cout << std::endl;
556  // throw std::invalid_argument("not enough test runs");
557
558  return sum/n - n;
559}
560template<class RandomAccessContainer>
561class generic_counter
562{
563public:
564  explicit generic_counter(unsigned int classes) : container(classes, 0) { }
565  void operator()(int i)
566  {
567    assert(i >= 0);
568    assert(static_cast<unsigned int>(i) < container.size());
569    ++container[i];
570  }
571  typename RandomAccessContainer::const_iterator begin() const 
572  { return container.begin(); }
573  typename RandomAccessContainer::const_iterator end() const 
574  { return container.end(); }
575
576private:
577  RandomAccessContainer container;
578};
579
580// chi_square test
581template<class Experiment, class Generator>
582double run_experiment(const Experiment & experiment, Generator gen, int n)
583{
584  generic_counter<std::vector<int> > v(experiment.classes());
585  experiment.run(gen, v, n);
586  return chi_square_value(v.begin(), v.end(),
587                          std::bind1st(std::mem_fun_ref(&Experiment::probability), 
588                                       experiment));
589}
590
591// number generator with experiment results (for nesting)
592template<class Experiment, class Generator>
593class experiment_generator_t
594{
595public:
596  experiment_generator_t(const Experiment & exper, Generator & gen, int n)
597    : experiment(exper), generator(gen), n(n) { }
598  double operator()() { return run_experiment(experiment, generator, n); }
599private:
600  const Experiment & experiment;
601  Generator & generator;
602  int n;
603};
604
605template<class Experiment, class Generator>
606experiment_generator_t<Experiment, Generator>
607experiment_generator(const Experiment & e, Generator & gen, int n)
608{
609  return experiment_generator_t<Experiment, Generator>(e, gen, n);
610}
611
612
613template<class Experiment, class Generator, class Distribution>
614class ks_experiment_generator_t
615{
616public:
617  ks_experiment_generator_t(const Experiment & exper, Generator & gen,
618                            const Distribution & distrib)
619    : experiment(exper), generator(gen), distribution(distrib) { }
620  double operator()() { return experiment.run(generator, distribution); }
621private:
622  const Experiment & experiment;
623  Generator & generator;
624  Distribution distribution;
625};
626
627template<class Experiment, class Generator, class Distribution>
628ks_experiment_generator_t<Experiment, Generator, Distribution>
629ks_experiment_generator(const Experiment & e, Generator & gen,
630                        const Distribution & distrib)
631{
632  return ks_experiment_generator_t<Experiment, Generator, Distribution>
633    (e, gen, distrib);
634}
635
636
637#endif /* STATISTIC_TESTS_HPP */
638
Note: See TracBrowser for help on using the repository browser.