core: (fixes #540) Add Bernoulli and Binomial Distributions

- Add BernoulliRandomVariable and BinomialRandomVariable to the random-variable-stream model
- Add BernoulliRandomVariable and BinomialRandomVariable to random-variables.rst
- Add a stanza to main-random-variable-stream.cc for BernoulliRandomVariable and BinomialRandomVariable
- Add tests for BernoulliRandomVariable and BinomialRandomVariable in random-variable-stream-test-suite.cc
- Update RELEASE_NOTES.md and CHANGES.md
This commit is contained in:
AlessioBugetti
2024-01-16 19:48:13 +01:00
parent 85d03d2353
commit 567566c259
8 changed files with 663 additions and 0 deletions

View File

@@ -333,6 +333,7 @@ Theodore Zhang (harryzwh@gmail.com)
Dizhi Zhou (dizhi.zhou@gmail.com)
Tolik Zinovyev (tolik@bu.edu)
Tommaso Zugno (tommasozugno@gmail.com)
Alessio Bugetti (alessiobugetti98@gmail.com)
hax0kartik (GCI 2019)
howie (GCI 2019)
InquisitivePenguin (GCI 2019)

View File

@@ -20,6 +20,7 @@ Changes from ns-3.40 to ns-3-dev
* (spectrum) `SpectrumSignalParameters` is extended to include two new members called: `spectrumChannelMatrix` and `precodingMatrix` which are the key information needed to support MIMO simulations.
* (wifi) Added new attribute `ChannelAccessManager:GenerateBackoffIfTxopWithoutTx` to invoke the backoff procedure when an AC gains the right to start a TXOP but it does not transmit any frame, provided that the queue is not actually empty. No transmission may occur,e.g., due to constraints associated with EMLSR operations. This possibility is specified by the current draft revision of the IEEE 802.11 standard.
* (core) Added `BernoulliRandomVariable` class implementing the bernoulli random variable, and `BinomialRandomVariable` class implementing the binomial random variable.
### Changes to existing API

View File

@@ -30,6 +30,7 @@ Release 3-dev
- (wifi) - Added EHT support for Ideal rate manager
- (wifi) - Reduce error rate model precision to fix infinite loop when Ideal rate manager is used with EHT
- (lr-wpan) !1794 - Group MAC primitives status enumerations into a single enumeration
- (core) !1802 - Added support for Bernoulli and Binomial random variables (`BernoulliRandomVariable`, `BinomialRandomVariable`)
### Bugs fixed

View File

@@ -261,6 +261,8 @@ can also create their own custom random variables by deriving from class
* class :cpp:class:`ZetaRandomVariable`
* class :cpp:class:`DeterministicRandomVariable`
* class :cpp:class:`EmpiricalRandomVariable`
* class :cpp:class:`BinomialRandomVariable`
* class :cpp:class:`BernoulliRandomVariable`
Semantics of RandomVariableStream objects
*****************************************

View File

@@ -555,6 +555,37 @@ main(int argc, char* argv[])
std::cout << "done" << std::endl;
}
{
std::cout << "BinomialRandomVariable......." << std::flush;
Gnuplot plot;
plot.SetTitle("BinomialRandomVariable");
plot.AppendExtra("set yrange [0:10]");
auto x = CreateObject<BinomialRandomVariable>();
x->SetAttribute("Trials", IntegerValue(10));
x->SetAttribute("Probability", DoubleValue(0.5));
plot.AddDataset(Histogram(x, probes, precision, "BinomialRandomVariable n=10 p=0.5"));
gnuplots.AddPlot(plot);
std::cout << "done" << std::endl;
}
{
std::cout << "BernoulliRandomVariable......." << std::flush;
Gnuplot plot;
plot.SetTitle("BernoulliRandomVariable");
plot.AppendExtra("set yrange [0:1]");
auto x = CreateObject<BernoulliRandomVariable>();
x->SetAttribute("Probability", DoubleValue(0.5));
plot.AddDataset(Histogram(x, probes, precision, "BernoulliRandomVariable p=0.5"));
gnuplots.AddPlot(plot);
std::cout << "done" << std::endl;
}
{
std::string gnuFile = plotFile + ".plt";
std::cout << "Writing Gnuplot file: " << gnuFile << "..." << std::flush;

View File

@@ -1738,4 +1738,123 @@ EmpiricalRandomVariable::Validate()
m_validated = true;
}
NS_OBJECT_ENSURE_REGISTERED(BinomialRandomVariable);
TypeId
BinomialRandomVariable::GetTypeId()
{
static TypeId tid =
TypeId("ns3::BinomialRandomVariable")
.SetParent<RandomVariableStream>()
.SetGroupName("Core")
.AddConstructor<BinomialRandomVariable>()
.AddAttribute("Trials",
"The number of trials.",
IntegerValue(10),
MakeIntegerAccessor(&BinomialRandomVariable::m_trials),
MakeIntegerChecker<uint32_t>(0))
.AddAttribute("Probability",
"The probability of success in each trial.",
DoubleValue(0.5),
MakeDoubleAccessor(&BinomialRandomVariable::m_probability),
MakeDoubleChecker<double>(0));
return tid;
}
BinomialRandomVariable::BinomialRandomVariable()
{
// m_trials and m_probability are initialized after constructor by attributes
NS_LOG_FUNCTION(this);
}
double
BinomialRandomVariable::GetValue(uint32_t trials, double probability)
{
NS_LOG_FUNCTION(this << trials << probability);
double successes = 0;
for (uint32_t i = 0; i < trials; ++i)
{
double v = Peek()->RandU01();
if (IsAntithetic())
{
v = (1 - v);
}
if (v <= probability)
{
successes += 1;
}
}
return successes;
}
uint32_t
BinomialRandomVariable::GetInteger(uint32_t trials, uint32_t probability)
{
NS_LOG_FUNCTION(this << trials << probability);
return static_cast<uint32_t>(GetValue(trials, probability));
}
double
BinomialRandomVariable::GetValue()
{
NS_LOG_FUNCTION(this);
return GetValue(m_trials, m_probability);
}
NS_OBJECT_ENSURE_REGISTERED(BernoulliRandomVariable);
TypeId
BernoulliRandomVariable::GetTypeId()
{
static TypeId tid =
TypeId("ns3::BernoulliRandomVariable")
.SetParent<RandomVariableStream>()
.SetGroupName("Core")
.AddConstructor<BernoulliRandomVariable>()
.AddAttribute("Probability",
"The probability of the random variable returning a value of 1.",
DoubleValue(0.5),
MakeDoubleAccessor(&BernoulliRandomVariable::m_probability),
MakeDoubleChecker<double>(0));
return tid;
}
BernoulliRandomVariable::BernoulliRandomVariable()
{
// m_probability is initialized after constructor by attributes
NS_LOG_FUNCTION(this);
}
double
BernoulliRandomVariable::GetValue(double probability)
{
NS_LOG_FUNCTION(this << probability);
double v = Peek()->RandU01();
if (IsAntithetic())
{
v = (1 - v);
}
return (v <= probability) ? 1.0 : 0.0;
}
uint32_t
BernoulliRandomVariable::GetInteger(uint32_t probability)
{
NS_LOG_FUNCTION(this << probability);
return static_cast<uint32_t>(GetValue(probability));
}
double
BernoulliRandomVariable::GetValue()
{
NS_LOG_FUNCTION(this);
return GetValue(m_probability);
}
} // namespace ns3

View File

@@ -18,6 +18,7 @@
* Authors: Rajib Bhattacharjea<raj.b@gatech.edu>
* Hadi Arbabi<marbabi@cs.odu.edu>
* Mathieu Lacage <mathieu.lacage@gmail.com>
* Alessio Bugetti <alessiobugetti98@gmail.com>
*
* Modified by Mitch Watrous <watrous@u.washington.edu>
*
@@ -2074,6 +2075,202 @@ class EmpiricalRandomVariable : public RandomVariableStream
}; // class EmpiricalRandomVariable
/**
* \ingroup randomvariable
* \brief The binomial distribution Random Number Generator (RNG).
*
* This class supports the creation of objects that return random numbers
* from a fixed binomial distribution. It also supports the generation of
* single random numbers from various binomial distributions.
*
* The probability mass function of a binomial variable
* is defined as:
*
* \f[
* P(k; n, p) = \binom{n}{k} p^k (1-p)^{n-k}, \\
* \quad k \in [0, n]
* \f]
*
* where \f$ n \f$ is the number of trials and \f$ p \f$ is the probability
* of success in each trial. The mean of this distribution
* is \f$ \mu = np \f$ and the variance is \f$ \sigma^2 = np(1-p) \f$.
*
* The Binomial RNG value \f$n\f$ for a given number of trials \f$ n \f$
* and success probability \f$ p \f$ is generated by
*
* \f[
* k = \sum_{i=1}^{n} I(u_i \leq p)
* \f]
*
* where \f$u_i\f$ is a uniform random variable on [0,1) for each trial,
* and \f$I\f$ is an indicator function that is 1 if \f$u_i \leq p\f$ and 0 otherwise.
* The sum of these indicator functions over all trials gives the total
* number of successes, which is the value of the binomial random variable.
*
* \par Example
*
* Here is an example of how to use this class:
* \code{.cc}
* uint32_t trials = 10;
* double probability = 0.5;
*
* Ptr<BinomialRandomVariable> x = CreateObject<BinomialRandomVariable> ();
* x->SetAttribute ("Trials", UintegerValue (trials));
* x->SetAttribute ("Probability", DoubleValue (probability));
*
* double successes = x->GetValue ();
* \endcode
*
* \par Antithetic Values.
*
* If an instance of this RNG is configured to return antithetic values,
* the actual value returned, \f$n'\f$, for the Binomial process is determined by:
*
* \f[
* k' = \sum_{i=1}^{n} I((1 - u_i) \leq p)
* \f]
*
* where \f$u_i\f$ is a uniform random variable on [0,1) for each trial.
* The antithetic approach uses \f$(1 - u_i)\f$ instead of \f$u_i\f$ in the indicator function.
*
* @par Efficiency and Alternative Methods
*
* There are alternative methods for generating binomial distributions that may offer greater
* efficiency. However, this implementation opts for a simpler approach. Although not as efficient,
* this method was chosen for its simplicity and sufficiency in most applications.
*/
class BinomialRandomVariable : public RandomVariableStream
{
public:
/**
* \brief Register this type.
* \return The object TypeId.
*/
static TypeId GetTypeId();
BinomialRandomVariable();
/**
* \copydoc GetValue()
* \param [in] trials Number of trials.
* \param [in] probability Probability of success in each trial.
* \return Returns a number within the range [0, trials] indicating the number of successful
* trials.
*/
double GetValue(uint32_t trials, double probability);
/**
* \copydoc GetValue(uint32_t,double)
* This function is similar to GetValue(), but it returns a uint32_t instead of a double.
*/
uint32_t GetInteger(uint32_t trials, uint32_t probability);
// Inherited
double GetValue() override;
using RandomVariableStream::GetInteger;
private:
/** The number of trials. */
uint32_t m_trials;
/** The probability of success in each trial. */
double m_probability;
}; // class BinomialRandomVariable
/**
* \ingroup randomvariable
* \brief The Bernoulli distribution Random Number Generator (RNG).
*
* This class supports the creation of objects that return random numbers
* from a fixed Bernoulli distribution. It also supports the generation of
* single random numbers from various Bernoulli distributions.
*
* The probability mass function of a Bernoulli variable
* is defined as:
*
* \f[
* P(n; p) = p^n (1-p)^{1-n}, \\
* \quad n \in \{0, 1\}
* \f]
*
* where \f$ p \f$ is the probability of success and n is the indicator of success, with n=1
* representing success and n=0 representing failure. The mean of this distribution is \f$ \mu = p
* \f$ and the variance is \f$ \sigma^2 = p(1-p) \f$.
*
* The Bernoulli RNG value \f$n\f$ is generated by
*
* \f[
* n =
* \begin{cases}
* 1 & \text{if } u \leq p \\
* 0 & \text{otherwise}
* \end{cases}
* \f]
*
* where \f$u\f$ is a uniform random variable on [0,1) and \f$p\f$ is the probability of success.
*
* \par Example
*
* Here is an example of how to use this class:
* \code{.cc}
* double probability = 0.5;
*
* Ptr<BernoulliRandomVariable> x = CreateObject<BernoulliRandomVariable> ();
* x->SetAttribute ("Probability", DoubleValue (probability));
*
* double success = x->GetValue ();
* \endcode
*
* \par Antithetic Values.
*
* If an instance of this RNG is configured to return antithetic values,
* the actual value returned, \f$x'\f$, is determined by:
*
* \f[
* x' =
* \begin{cases}
* 1 & \text{if } (1 - u) \leq p \\
* 0 & \text{otherwise}
* \end{cases}
* \f]
*
* where \f$u\f$ is a uniform random variable on [0,1) and \f$p\f$ is the probability of success.
*/
class BernoulliRandomVariable : public RandomVariableStream
{
public:
/**
* \brief Register this type.
* \return The object TypeId.
*/
static TypeId GetTypeId();
BernoulliRandomVariable();
/**
* \copydoc GetValue()
* \param [in] probability Probability of success.
* \return Returns 1 if the trial is successful, or 0 if it is not.
*/
double GetValue(double probability);
/**
* \copydoc GetValue(double)
* This function is similar to GetValue(), but it returns a uint32_t instead of a double.
*/
uint32_t GetInteger(uint32_t probability);
// Inherited
double GetValue() override;
using RandomVariableStream::GetInteger;
private:
/** The probability of success. */
double m_probability;
}; // class BernoulliRandomVariable
} // namespace ns3
#endif /* RANDOM_VARIABLE_STREAM_H */

View File

@@ -34,6 +34,7 @@
#include <fstream>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_zeta.h>
using namespace ns3;
@@ -2565,6 +2566,312 @@ NormalCachingTestCase::DoRun()
NS_TEST_ASSERT_MSG_GT(v2, 0, "Incorrect value returned, expected > 0");
}
/**
* \ingroup rng-tests
* Test case for bernoulli distribution random variable stream generator
*/
class BernoulliTestCase : public TestCaseBase
{
public:
// Constructor
BernoulliTestCase();
// Inherited
double ChiSquaredTest(Ptr<RandomVariableStream> rng) const override;
private:
// Inherited
void DoRun() override;
/** Tolerance for testing rng values against expectation, in rms. */
static constexpr double TOLERANCE{5};
};
BernoulliTestCase::BernoulliTestCase()
: TestCaseBase("Bernoulli Random Variable Stream Generator")
{
}
double
BernoulliTestCase::ChiSquaredTest(Ptr<RandomVariableStream> rng) const
{
gsl_histogram* h = gsl_histogram_alloc(2);
auto range = UniformHistogramBins(h, 0, 1);
double p = 0.5;
std::vector<double> expected = {N_MEASUREMENTS * (1 - p), N_MEASUREMENTS * p};
double chiSquared = ChiSquared(h, expected, rng);
gsl_histogram_free(h);
return chiSquared;
}
void
BernoulliTestCase::DoRun()
{
NS_LOG_FUNCTION(this);
SetTestSuiteSeed();
auto generator = RngGenerator<BernoulliRandomVariable>();
double sum = ChiSquaredsAverage(&generator, N_RUNS);
double maxStatistic = gsl_cdf_chisq_Qinv(0.05, 1);
NS_TEST_ASSERT_MSG_LT(sum, maxStatistic, "Chi-squared statistic out of range");
double probability = 0.5;
// Create the RNG with the specified range.
Ptr<BernoulliRandomVariable> x = CreateObject<BernoulliRandomVariable>();
x->SetAttribute("Probability", DoubleValue(probability));
// Calculate the mean of these values.
double mean = probability;
double valueMean = Average(x);
double expectedMean = mean;
double expectedRms = std::sqrt(mean / N_MEASUREMENTS);
// Test that values have approximately the right mean value.
NS_TEST_ASSERT_MSG_EQ_TOL(valueMean,
expectedMean,
expectedRms * TOLERANCE,
"Wrong mean value.");
}
/**
* \ingroup rng-tests
* Test case for antithetic bernoulli distribution random variable stream generator
*/
class BernoulliAntitheticTestCase : public TestCaseBase
{
public:
// Constructor
BernoulliAntitheticTestCase();
// Inherited
double ChiSquaredTest(Ptr<RandomVariableStream> rng) const override;
private:
// Inherited
void DoRun() override;
/** Tolerance for testing rng values against expectation, in rms. */
static constexpr double TOLERANCE{5};
};
BernoulliAntitheticTestCase::BernoulliAntitheticTestCase()
: TestCaseBase("Antithetic Bernoulli Random Variable Stream Generator")
{
}
double
BernoulliAntitheticTestCase::ChiSquaredTest(Ptr<RandomVariableStream> rng) const
{
gsl_histogram* h = gsl_histogram_alloc(2);
auto range = UniformHistogramBins(h, 0, 1);
double p = 0.5;
std::vector<double> expected = {N_MEASUREMENTS * (1 - p), N_MEASUREMENTS * p};
double chiSquared = ChiSquared(h, expected, rng);
gsl_histogram_free(h);
return chiSquared;
}
void
BernoulliAntitheticTestCase::DoRun()
{
NS_LOG_FUNCTION(this);
SetTestSuiteSeed();
auto generator = RngGenerator<BernoulliRandomVariable>(true);
double sum = ChiSquaredsAverage(&generator, N_RUNS);
double maxStatistic = gsl_cdf_chisq_Qinv(0.05, 1);
NS_TEST_ASSERT_MSG_LT(sum, maxStatistic, "Chi-squared statistic out of range");
double probability = 0.5;
// Create the RNG with the specified range.
Ptr<BernoulliRandomVariable> x = CreateObject<BernoulliRandomVariable>();
x->SetAttribute("Probability", DoubleValue(probability));
// Make this generate antithetic values.
x->SetAttribute("Antithetic", BooleanValue(true));
// Calculate the mean of these values.
double mean = probability;
double valueMean = Average(x);
double expectedMean = mean;
double expectedRms = std::sqrt(mean / N_MEASUREMENTS);
// Test that values have approximately the right mean value.
NS_TEST_ASSERT_MSG_EQ_TOL(valueMean,
expectedMean,
expectedRms * TOLERANCE,
"Wrong mean value.");
}
/**
* \ingroup rng-tests
* Test case for binomial distribution random variable stream generator
*/
class BinomialTestCase : public TestCaseBase
{
public:
// Constructor
BinomialTestCase();
// Inherited
double ChiSquaredTest(Ptr<RandomVariableStream> rng) const override;
private:
// Inherited
void DoRun() override;
/** Tolerance for testing rng values against expectation, in rms. */
static constexpr double TOLERANCE{5};
};
BinomialTestCase::BinomialTestCase()
: TestCaseBase("Binomial Random Variable Stream Generator")
{
}
double
BinomialTestCase::ChiSquaredTest(Ptr<RandomVariableStream> rng) const
{
uint32_t trials = 10;
double probability = 0.5;
gsl_histogram* h = gsl_histogram_alloc(trials + 1);
auto range = UniformHistogramBins(h, 0, trials);
std::vector<double> expected(trials + 1);
for (std::size_t i = 0; i < trials + 1; ++i)
{
expected[i] = N_MEASUREMENTS * gsl_ran_binomial_pdf(i, probability, trials);
}
double chiSquared = ChiSquared(h, expected, rng);
gsl_histogram_free(h);
return chiSquared;
}
void
BinomialTestCase::DoRun()
{
NS_LOG_FUNCTION(this);
SetTestSuiteSeed();
uint32_t trials = 10;
double probability = 0.5;
auto generator = RngGenerator<BinomialRandomVariable>();
double sum = ChiSquaredsAverage(&generator, N_RUNS);
double maxStatistic = gsl_cdf_chisq_Qinv(0.05, trials);
NS_TEST_ASSERT_MSG_LT(sum, maxStatistic, "Chi-squared statistic out of range");
// Create the RNG with the specified range.
Ptr<BinomialRandomVariable> x = CreateObject<BinomialRandomVariable>();
x->SetAttribute("Trials", IntegerValue(trials));
x->SetAttribute("Probability", DoubleValue(probability));
// Calculate the mean of these values.
double mean = trials * probability;
double valueMean = Average(x);
double expectedMean = mean;
double expectedRms = std::sqrt(mean / N_MEASUREMENTS);
// Test that values have approximately the right mean value.
NS_TEST_ASSERT_MSG_EQ_TOL(valueMean,
expectedMean,
expectedRms * TOLERANCE,
"Wrong mean value.");
}
/**
* \ingroup rng-tests
* Test case for antithetic binomial distribution random variable stream generator
*/
class BinomialAntitheticTestCase : public TestCaseBase
{
public:
// Constructor
BinomialAntitheticTestCase();
// Inherited
double ChiSquaredTest(Ptr<RandomVariableStream> rng) const override;
private:
// Inherited
void DoRun() override;
/** Tolerance for testing rng values against expectation, in rms. */
static constexpr double TOLERANCE{5};
};
BinomialAntitheticTestCase::BinomialAntitheticTestCase()
: TestCaseBase("Antithetic Binomial Random Variable Stream Generator")
{
}
double
BinomialAntitheticTestCase::ChiSquaredTest(Ptr<RandomVariableStream> rng) const
{
uint32_t trials = 10;
double probability = 0.5;
gsl_histogram* h = gsl_histogram_alloc(trials + 1);
auto range = UniformHistogramBins(h, 0, trials);
std::vector<double> expected(trials + 1);
for (std::size_t i = 0; i < trials + 1; ++i)
{
expected[i] = N_MEASUREMENTS * gsl_ran_binomial_pdf(i, probability, trials);
}
double chiSquared = ChiSquared(h, expected, rng);
gsl_histogram_free(h);
return chiSquared;
}
void
BinomialAntitheticTestCase::DoRun()
{
NS_LOG_FUNCTION(this);
SetTestSuiteSeed();
uint32_t trials = 10;
double probability = 0.5;
auto generator = RngGenerator<BinomialRandomVariable>(true);
double sum = ChiSquaredsAverage(&generator, N_RUNS);
double maxStatistic = gsl_cdf_chisq_Qinv(0.05, trials);
NS_TEST_ASSERT_MSG_LT(sum, maxStatistic, "Chi-squared statistic out of range");
// Create the RNG with the specified range.
Ptr<BinomialRandomVariable> x = CreateObject<BinomialRandomVariable>();
x->SetAttribute("Trials", IntegerValue(trials));
x->SetAttribute("Probability", DoubleValue(probability));
// Make this generate antithetic values.
x->SetAttribute("Antithetic", BooleanValue(true));
// Calculate the mean of these values.
double mean = trials * probability;
double valueMean = Average(x);
double expectedMean = mean;
double expectedRms = std::sqrt(mean / N_MEASUREMENTS);
// Test that values have approximately the right mean value.
NS_TEST_ASSERT_MSG_EQ_TOL(valueMean,
expectedMean,
expectedRms * TOLERANCE,
"Wrong mean value.");
}
/**
* \ingroup rng-tests
* RandomVariableStream test suite, covering all random number variable
@@ -2617,6 +2924,10 @@ RandomVariableSuite::RandomVariableSuite()
AddTestCase(new EmpiricalAntitheticTestCase);
/// Issue #302: NormalRandomVariable produces stale values
AddTestCase(new NormalCachingTestCase);
AddTestCase(new BernoulliTestCase);
AddTestCase(new BernoulliAntitheticTestCase);
AddTestCase(new BinomialTestCase);
AddTestCase(new BinomialAntitheticTestCase);
}
static RandomVariableSuite randomVariableSuite; //!< Static variable for test initialization