feat: add geometric distribution (#1699)

* feat: add geometric distribution * test: Add test for geometric dist * fix: Make range_tries inclusive * docs: Add documentation for geometric distribution * fix: clang warnings in geometric_dist * updating DIRECTORY.md * fix: Remove extra line Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Remove file name Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Add return value of void function Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Add comment for test function Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Update successful test message Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Add geometric_dist namespace Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Remove extra line Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Close comment * docs: Fix documentation Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Add const to parameter Co-authored-by: David Leal <halfpacho@gmail.com> * fix: Make class methods const Co-authored-by: github-actions <${GITHUB_ACTOR}@users.noreply.github.com> Co-authored-by: David Leal <halfpacho@gmail.com>
2023-10-11 13:05:55 +08:00 · 2021-10-13 21:22:32 +00:00 · 2021-10-13 21:22:32 +00:00 · 5a654fb85b
commit 5a654fb85b
parent 89bfc1244e
2 changed files with 237 additions and 0 deletions
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@ -266,6 +266,7 @@
  * [Addition Rule](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/probability/addition_rule.cpp)
  * [Bayes Theorem](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/probability/bayes_theorem.cpp)
  * [Binomial Dist](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/probability/binomial_dist.cpp)
+  * [Geometric Dist](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/probability/geometric_dist.cpp)
  * [Poisson Dist](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/probability/poisson_dist.cpp)
  * [Windowed Median](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/probability/windowed_median.cpp)

--- a/probability/geometric_dist.cpp
+++ b/probability/geometric_dist.cpp
@ -0,0 +1,236 @@
+/**
+ * @file
+ * @brief [Geometric Distribution](https://en.wikipedia.org/wiki/Geometric_distribution)
+ *
+ * @details
+ * The geometric distribution models the experiment of doing Bernoulli trials until a
+ * sucess was observed. There are two formulations of the geometric distribution:
+ * 1) The probability distribution of the number X of Bernoulli trials needed to get one success, supported on the set { 1, 2, 3, ... }
+ * 2) The probability distribution of the number Y = X − 1 of failures before the first success, supported on the set { 0, 1, 2, 3, ... }
+ * Here, the first one is implemented.
+ *
+ * Common variables used:
+ * p - The success probability
+ * k - The number of tries
+ *
+ * @author [Domenic Zingsheim](https://github.com/DerAndereDomenic)
+ */
+
+#include <cassert>    /// for assert
+#include <cmath>      /// for math functions
+#include <cstdint>    /// for fixed size data types
+#include <ctime>      /// for time to initialize rng
+#include <iostream>   /// for std::cout
+#include <limits>     /// for std::numeric_limits
+#include <random>     /// for random numbers
+#include <vector>     /// for std::vector
+
+/**
+ * @namespace probability
+ * @brief Probability algorithms
+ */
+namespace probability {
+/**
+ * @namespace geometric_dist
+ * @brief Functions for the [Geometric Distribution](https://en.wikipedia.org/wiki/Geometric_distribution) algorithm implementation
+ */
+namespace geometric_dist {
+/**
+ * @brief Returns a random number between [0,1]
+ * @returns A uniformly distributed random number between 0 (included) and 1 (included)
+ */
+float generate_uniform() {
+    return static_cast<float>(rand()) / static_cast<float>(RAND_MAX);
+}
+
+/**
+ * @brief A class to model the geometric distribution
+ */
+class geometric_distribution
+{
+private:
+    float p;    ///< The succes probability p
+
+public:
+    /**
+     * @brief Constructor for the geometric distribution
+     * @param p The success probability
+     */
+    explicit geometric_distribution(const float& p) : p(p) {}
+
+    /**
+     * @brief The expected value of a geometrically distributed random variable X
+     * @returns E[X] = 1/p
+     */
+    float expected_value() const {
+       return 1.0f/ p;
+    }
+
+    /**
+     * @brief The variance of a geometrically distributed random variable X
+     * @returns V[X] = (1 - p) / p^2
+     */
+    float variance() const {
+       return (1.0f - p) / (p * p);
+    }
+
+    /**
+     * @brief The standard deviation of a geometrically distributed random variable X
+     * @returns \sigma = \sqrt{V[X]}
+     */
+    float standard_deviation() const {
+       return std::sqrt(variance());
+    }
+
+    /**
+     * @brief The probability density function
+     * @details As we use the first definition of the geometric series (1),
+     * we are doing k - 1 failed trials and the k-th trial is a success.
+     * @param k The number of trials to observe the first success in [1,\infty)
+     * @returns A number between [0,1] according to p * (1-p)^{k-1}
+     */
+    float probability_density(const uint32_t& k) const {
+        return std::pow((1.0f - p), static_cast<float>(k - 1)) * p;
+    }
+
+    /**
+     * @brief The cumulative distribution function
+     * @details The sum of all probabilities up to (and including) k trials. Basically CDF(k) = P(x <= k)
+     * @param k The number of trials in [1,\infty)
+     * @returns The probability to have success within k trials
+     */
+    float cumulative_distribution(const uint32_t& k) const {
+        return 1.0f - std::pow((1.0f - p), static_cast<float>(k));
+    }
+
+    /**
+     * @brief The inverse cumulative distribution function
+     * @details This functions answers the question: Up to how many trials are needed to have success with a probability of cdf?
+     * The exact floating point value is reported.
+     * @param cdf The probability in [0,1]
+     * @returns The number of (exact) trials.
+     */
+    float inverse_cumulative_distribution(const float& cdf) const {
+        return std::log(1.0f - cdf) / std::log(1.0f - p);
+    }
+
+    /**
+     * @brief Generates a (discrete) sample according to the geometrical distribution
+     * @returns A geometrically distributed number in [1,\infty)
+     */
+    uint32_t draw_sample() const {
+        float uniform_sample = generate_uniform();
+        return static_cast<uint32_t>(inverse_cumulative_distribution(uniform_sample)) + 1;
+    }
+
+    /**
+     * @brief This function computes the probability to have success in a given range of tries
+     * @details Computes P(min_tries <= x <= max_tries).
+     * Can be used to calculate P(x >= min_tries) by not passing a second argument.
+     * Can be used to calculate P(x <= max_tries) by passing 1 as the first argument
+     * @param min_tries The minimum number of tries in [1,\infty) (inclusive)
+     * @param max_tries The maximum number of tries in [min_tries, \infty) (inclusive)
+     * @returns The probability of having success within a range of tries [min_tries, max_tries]
+     */
+    float range_tries(const uint32_t& min_tries = 1, const uint32_t& max_tries = std::numeric_limits<uint32_t>::max()) const {
+        float cdf_lower = cumulative_distribution(min_tries - 1);
+        float cdf_upper = max_tries == std::numeric_limits<uint32_t>::max() ? 1.0f : cumulative_distribution(max_tries);
+        return cdf_upper - cdf_lower;
+    }
+};
+}  // namespace geometric_dist
+}  // namespace probability
+
+/**
+ * @brief Tests the sampling method of the geometric distribution
+ * @details Draws 1000000 random samples and estimates mean and variance
+ * These should be close to the expected value and variance of the given distribution to pass.
+ * @param dist The distribution to test
+ */
+void sample_test(const probability::geometric_dist::geometric_distribution& dist) {
+    uint32_t n_tries = 1000000;
+    std::vector<float> tries;
+    tries.resize(n_tries);
+
+    float mean = 0.0f;
+    for (uint32_t i = 0; i < n_tries; ++i) {
+        tries[i] = static_cast<float>(dist.draw_sample());
+        mean += tries[i];
+    }
+
+    mean /= static_cast<float>(n_tries);
+
+    float var = 0.0f;
+    for (uint32_t i = 0; i < n_tries; ++i) {
+        var += (tries[i] - mean) * (tries[i] - mean);
+    }
+
+    //Unbiased estimate of variance
+    var /= static_cast<float>(n_tries - 1);
+
+    std::cout << "This value should be near " << dist.expected_value() << ": " << mean << std::endl;
+    std::cout << "This value should be near " << dist.variance() << ": " << var << std::endl;
+}
+
+/**
+ * @brief Self-test implementations
+ * @returns void
+ */
+static void test() {
+    probability::geometric_dist::geometric_distribution dist(0.3);
+
+    const float threshold = 1e-3f;
+
+    std::cout << "Starting tests for p = 0.3..." << std::endl;
+    assert(std::abs(dist.expected_value() - 3.33333333f) < threshold);
+    assert(std::abs(dist.variance() - 7.77777777f) < threshold);
+    assert(std::abs(dist.standard_deviation() - 2.788866755) < threshold);
+    assert(std::abs(dist.probability_density(5) - 0.07203) < threshold);
+    assert(std::abs(dist.cumulative_distribution(6) - 0.882351) < threshold);
+    assert(std::abs(dist.inverse_cumulative_distribution(dist.cumulative_distribution(8)) - 8) < threshold);
+    assert(std::abs(dist.range_tries() - 1.0f) < threshold);
+    assert(std::abs(dist.range_tries(3) - 0.49f) < threshold);
+    assert(std::abs(dist.range_tries(5, 11) - 0.2203267f) < threshold);
+    std::cout << "All tests passed" << std::endl;
+    sample_test(dist);
+
+    dist = probability::geometric_dist::geometric_distribution(0.5f);
+
+    std::cout << "Starting tests for p = 0.5..." << std::endl;
+    assert(std::abs(dist.expected_value() - 2.0f) < threshold);
+    assert(std::abs(dist.variance() - 2.0f) < threshold);
+    assert(std::abs(dist.standard_deviation() - 1.4142135f) < threshold);
+    assert(std::abs(dist.probability_density(5) - 0.03125) < threshold);
+    assert(std::abs(dist.cumulative_distribution(6) - 0.984375) < threshold);
+    assert(std::abs(dist.inverse_cumulative_distribution(dist.cumulative_distribution(8)) - 8) < threshold);
+    assert(std::abs(dist.range_tries() - 1.0f) < threshold);
+    assert(std::abs(dist.range_tries(3) - 0.25f) < threshold);
+    assert(std::abs(dist.range_tries(5, 11) - 0.062011f) < threshold);
+    std::cout << "All tests passed" << std::endl;
+    sample_test(dist);
+
+    dist = probability::geometric_dist::geometric_distribution(0.8f);
+
+    std::cout << "Starting tests for p = 0.8..." << std::endl;
+    assert(std::abs(dist.expected_value() - 1.25f) < threshold);
+    assert(std::abs(dist.variance() - 0.3125f) < threshold);
+    assert(std::abs(dist.standard_deviation() - 0.559016f) < threshold);
+    assert(std::abs(dist.probability_density(5) - 0.00128) < threshold);
+    assert(std::abs(dist.cumulative_distribution(6) - 0.999936) < threshold);
+    assert(std::abs(dist.inverse_cumulative_distribution(dist.cumulative_distribution(8)) - 8) < threshold);
+    assert(std::abs(dist.range_tries() - 1.0f) < threshold);
+    assert(std::abs(dist.range_tries(3) - 0.04f) < threshold);
+    assert(std::abs(dist.range_tries(5, 11) - 0.00159997f) < threshold);
+    std::cout << "All tests have successfully passed!" << std::endl;
+    sample_test(dist);
+}
+
+/**
+ * @brief Main function
+ * @return 0 on exit
+ */
+int main() {
+    srand(time(nullptr));
+    test();  // run self-test implementations
+    return 0;
+}