113b21bdaf
2011-03-24 Paolo Carlini <paolo.carlini@oracle.com> * include/bits/random.h (negative_binomial_distribution<>:: negative_binomial_distribution(_IntType, double), negative_binomial_distribution<>:: negative_binomial_distribution(const param_type&)): Fix construction of _M_gd. * include/bits/random.tcc (negative_binomial_distribution<>:: operator()): Fix computation, per Leger's algorithm. * testsuite/util/testsuite_random.h (discrete_pdf, negative_binomial_pdf, poisson_pdf, uniform_int_pdf): New. (binomial_pdf): Swap last two parameters. * testsuite/26_numerics/random/discrete_distribution/ operators/values.cc: New. * testsuite/26_numerics/random/negative_binomial_distribution/ operators/values.cc: Likewise. * testsuite/26_numerics/random/poisson_distribution/ operators/values.cc: Likewise. * testsuite/26_numerics/random/uniform_int_distribution/ operators/values.cc: Likewise. * testsuite/26_numerics/random/binomial_distribution/ operators/values.cc: Adjust. From-SVN: r171411
182 lines
3.7 KiB
C++
182 lines
3.7 KiB
C++
// -*- C++ -*-
|
|
|
|
// Copyright (C) 2011 Free Software Foundation, Inc.
|
|
//
|
|
// This file is part of the GNU ISO C++ Library. This library is free
|
|
// software; you can redistribute it and/or modify it under the terms
|
|
// of the GNU General Public License as published by the Free Software
|
|
// Foundation; either version 3, or (at your option) any later
|
|
// version.
|
|
|
|
// This library is distributed in the hope that it will be useful, but
|
|
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
// General Public License for more details.
|
|
|
|
// You should have received a copy of the GNU General Public License along
|
|
// with this library; see the file COPYING3. If not see
|
|
// <http://www.gnu.org/licenses/>.
|
|
|
|
/**
|
|
* @file testsuite_random.h
|
|
*/
|
|
|
|
#ifndef _GLIBCXX_TESTSUITE_RANDOM_H
|
|
#define _GLIBCXX_TESTSUITE_RANDOM_H
|
|
|
|
#include <cmath>
|
|
#include <initializer_list>
|
|
#include <testsuite_hooks.h>
|
|
|
|
namespace __gnu_test
|
|
{
|
|
// Adapted for libstdc++ from GNU gsl-1.14/randist/test.c
|
|
// Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2010
|
|
// James Theiler, Brian Gough
|
|
template<unsigned long BINS = 100,
|
|
unsigned long N = 100000,
|
|
typename Distribution, typename Pdf>
|
|
void
|
|
testDiscreteDist(Distribution& f, Pdf pdf)
|
|
{
|
|
bool test __attribute__((unused)) = true;
|
|
double count[BINS], p[BINS];
|
|
|
|
for (unsigned long i = 0; i < BINS; i++)
|
|
count[i] = 0;
|
|
|
|
for (unsigned long i = 0; i < N; i++)
|
|
{
|
|
auto r = f();
|
|
if (r >= 0 && r < BINS)
|
|
count[r]++;
|
|
}
|
|
|
|
for (unsigned long i = 0; i < BINS; i++)
|
|
p[i] = pdf(i);
|
|
|
|
for (unsigned long i = 0; i < BINS; i++)
|
|
{
|
|
bool status_i;
|
|
double d = std::abs(count[i] - N * p[i]);
|
|
|
|
if (p[i] != 0)
|
|
{
|
|
double s = d / std::sqrt(N * p[i]);
|
|
status_i = (s > 5) && (d > 1);
|
|
}
|
|
else
|
|
status_i = (count[i] != 0);
|
|
|
|
VERIFY( !status_i );
|
|
}
|
|
}
|
|
|
|
inline double
|
|
bernoulli_pdf(int k, double p)
|
|
{
|
|
if (k == 0)
|
|
return 1 - p;
|
|
else if (k == 1)
|
|
return p;
|
|
else
|
|
return 0.0;
|
|
}
|
|
|
|
#ifdef _GLIBCXX_USE_C99_MATH_TR1
|
|
inline double
|
|
binomial_pdf(int k, int n, double p)
|
|
{
|
|
if (k < 0 || k > n)
|
|
return 0.0;
|
|
else
|
|
{
|
|
double q;
|
|
|
|
if (p == 0.0)
|
|
q = (k == 0) ? 1.0 : 0.0;
|
|
else if (p == 1.0)
|
|
q = (k == n) ? 1.0 : 0.0;
|
|
else
|
|
{
|
|
double ln_Cnk = (std::lgamma(n + 1.0) - std::lgamma(k + 1.0)
|
|
- std::lgamma(n - k + 1.0));
|
|
q = ln_Cnk + k * std::log(p) + (n - k) * std::log1p(-p);
|
|
q = std::exp(q);
|
|
}
|
|
|
|
return q;
|
|
}
|
|
}
|
|
#endif
|
|
|
|
inline double
|
|
discrete_pdf(int k, std::initializer_list<double> wl)
|
|
{
|
|
if (!wl.size())
|
|
wl = { 1.0 };
|
|
|
|
if (k < 0 || k >= wl.size())
|
|
return 0.0;
|
|
else
|
|
{
|
|
double sum = 0.0;
|
|
for (auto it = wl.begin(); it != wl.end(); ++it)
|
|
sum += *it;
|
|
return wl.begin()[k] / sum;
|
|
}
|
|
}
|
|
|
|
inline double
|
|
geometric_pdf(int k, double p)
|
|
{
|
|
if (k < 0)
|
|
return 0.0;
|
|
else if (k == 0)
|
|
return p;
|
|
else
|
|
return p * std::pow(1 - p, k);
|
|
}
|
|
|
|
#ifdef _GLIBCXX_USE_C99_MATH_TR1
|
|
inline double
|
|
negative_binomial_pdf(int k, int n, double p)
|
|
{
|
|
if (k < 0)
|
|
return 0.0;
|
|
else
|
|
{
|
|
double f = std::lgamma(k + (double)n);
|
|
double a = std::lgamma(n);
|
|
double b = std::lgamma(k + 1.0);
|
|
|
|
return std::exp(f - a - b) * std::pow(p, n) * std::pow(1 - p, k);
|
|
}
|
|
}
|
|
|
|
inline double
|
|
poisson_pdf(int k, double mu)
|
|
{
|
|
if (k < 0)
|
|
return 0.0;
|
|
else
|
|
{
|
|
double lf = std::lgamma(k + 1.0);
|
|
return std::exp(std::log(mu) * k - lf - mu);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
inline double
|
|
uniform_int_pdf(int k, int a, int b)
|
|
{
|
|
if (k < 0 || k < a || k > b)
|
|
return 0.0;
|
|
else
|
|
return 1.0 / (b - a + 1.0);
|
|
}
|
|
|
|
} // namespace __gnu_test
|
|
|
|
#endif // #ifndef _GLIBCXX_TESTSUITE_RANDOM_H
|