diff --git a/AUTHORS b/AUTHORS index eb103c65e..6fa59a431 100644 --- a/AUTHORS +++ b/AUTHORS @@ -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) diff --git a/CHANGES.md b/CHANGES.md index 0759b60f1..964e7a31f 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -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 diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 9a869c65c..32562352b 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -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 diff --git a/doc/manual/source/random-variables.rst b/doc/manual/source/random-variables.rst index a42f9666a..ed6e9afc3 100644 --- a/doc/manual/source/random-variables.rst +++ b/doc/manual/source/random-variables.rst @@ -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 ***************************************** diff --git a/src/core/examples/main-random-variable-stream.cc b/src/core/examples/main-random-variable-stream.cc index 96ef3f8e5..ea435936b 100644 --- a/src/core/examples/main-random-variable-stream.cc +++ b/src/core/examples/main-random-variable-stream.cc @@ -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(); + 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(); + 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; diff --git a/src/core/model/random-variable-stream.cc b/src/core/model/random-variable-stream.cc index 229b96829..4364419ea 100644 --- a/src/core/model/random-variable-stream.cc +++ b/src/core/model/random-variable-stream.cc @@ -1738,4 +1738,123 @@ EmpiricalRandomVariable::Validate() m_validated = true; } +NS_OBJECT_ENSURE_REGISTERED(BinomialRandomVariable); + +TypeId +BinomialRandomVariable::GetTypeId() +{ + static TypeId tid = + TypeId("ns3::BinomialRandomVariable") + .SetParent() + .SetGroupName("Core") + .AddConstructor() + .AddAttribute("Trials", + "The number of trials.", + IntegerValue(10), + MakeIntegerAccessor(&BinomialRandomVariable::m_trials), + MakeIntegerChecker(0)) + .AddAttribute("Probability", + "The probability of success in each trial.", + DoubleValue(0.5), + MakeDoubleAccessor(&BinomialRandomVariable::m_probability), + MakeDoubleChecker(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(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() + .SetGroupName("Core") + .AddConstructor() + .AddAttribute("Probability", + "The probability of the random variable returning a value of 1.", + DoubleValue(0.5), + MakeDoubleAccessor(&BernoulliRandomVariable::m_probability), + MakeDoubleChecker(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(GetValue(probability)); +} + +double +BernoulliRandomVariable::GetValue() +{ + NS_LOG_FUNCTION(this); + return GetValue(m_probability); +} + } // namespace ns3 diff --git a/src/core/model/random-variable-stream.h b/src/core/model/random-variable-stream.h index e7b398995..7ddb6f9e7 100644 --- a/src/core/model/random-variable-stream.h +++ b/src/core/model/random-variable-stream.h @@ -18,6 +18,7 @@ * Authors: Rajib Bhattacharjea * Hadi Arbabi * Mathieu Lacage + * Alessio Bugetti * * Modified by Mitch Watrous * @@ -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 x = CreateObject (); + * 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 x = CreateObject (); + * 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 */ diff --git a/src/core/test/random-variable-stream-test-suite.cc b/src/core/test/random-variable-stream-test-suite.cc index d9e7c3d0a..f91e50a6d 100644 --- a/src/core/test/random-variable-stream-test-suite.cc +++ b/src/core/test/random-variable-stream-test-suite.cc @@ -34,6 +34,7 @@ #include #include #include +#include #include 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 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 rng) const +{ + gsl_histogram* h = gsl_histogram_alloc(2); + auto range = UniformHistogramBins(h, 0, 1); + + double p = 0.5; + std::vector 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(); + 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 x = CreateObject(); + 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 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 rng) const +{ + gsl_histogram* h = gsl_histogram_alloc(2); + auto range = UniformHistogramBins(h, 0, 1); + + double p = 0.5; + std::vector 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(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 x = CreateObject(); + 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 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 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 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(); + 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 x = CreateObject(); + 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 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 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 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(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 x = CreateObject(); + 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