refactor(binomial-distribution): replace batch sampler with Walker's alias table

Sources/StatKit/Descriptive Statistics/Distribution/Discrete/BinomialDistribution.swift

@@ -6,28 +6,28 @@ public struct BinomialDistribution: DiscreteDistribution, UnivariateDistribution
   ///
   /// This property models what is commonly referred to as `p`.
   public let probability: Double
-  
+
   /// The number of trials under study.
   ///
   /// This property models what is commonly referred to as `n`.
   public let trials: Int
-  
+
   public var mean: Double {
     Double(trials) * probability
   }
-  
+
   public var variance: Double {
     Double(trials) * probability * (1 - probability)
   }
-  
+
   public var skewness: Double {
     (1 - 2 * probability) / (variance).squareRoot()
   }
-  
+
   public var excessKurtosis: Double {
     (1 - 6 * probability * (1 - probability)) / variance
   }
-  
+
   /// Creates a Binomial Distribution for a specified number of trials and a probability.
   /// - parameter probability: The probability of a successful trial.
   /// - parameter trials: The number of trials under investigation.
@@ -42,81 +42,64 @@ public struct BinomialDistribution: DiscreteDistribution, UnivariateDistribution
       0 <= trials,
       "Possibilities need to non-negative (\(trials) provided)."
     )
-    
+
     self.probability = probability
     self.trials = trials
   }
-  
+
   public func pmf(x: Int, logarithmic: Bool = false) -> Double {
     let coefficient = Double(choose(n: trials, k: x))
     let successes: Double = .pow(probability, Double(x))
     let failures: Double = .pow(1 - probability, Double(trials - x))
     let result = coefficient * successes * failures
-    
+
     switch result {
       case ...0:
         return logarithmic ? -.infinity : 0
-      
+
       case 1...:
         return logarithmic ? 0 : 1
-        
+
       default:
         return logarithmic ? .log(result) : result
     }
   }
-  
+
   public func cdf(x: Int, logarithmic: Bool = false) -> Double {
     let result = (0 ... x).reduce(into: 0) { (result, successes) in
       result += pmf(x: successes)
     }
-    
+
     switch result {
       case ...0:
         return logarithmic ? -.infinity : 0
-      
+
       case 1...:
         return logarithmic ? 0 : 1
-        
+
       default:
         return logarithmic ? .log(result) : result
     }
   }
-  
+
   public func sample() -> Int {
     let probability = Double.random(in: 0 ... 1)
     var cumulativeProbability = 0.0
-    
+
     for successes in 0 ... trials {
       cumulativeProbability += pmf(x: successes)
       if probability <= cumulativeProbability {
         return successes
       }
     }
-    
+
     return trials
   }
-  
+
   public func sample(_ numberOfElements: Int) -> [Int] {
-    var uniformGenerator = Xoroshiro256StarStar()
-    let probabilities = (1 ... numberOfElements).map { _ in Double.random(in: 0 ... 1, using: &uniformGenerator) }
-    var cumulativeProbs = [Double](repeating: 0, count: trials + 1)
-    cumulativeProbs[0] = pmf(x: 0)
-
-    for x in 1 ... trials {
-      cumulativeProbs[x] = cumulativeProbs[x - 1] + pmf(x: x)
-    }
-
-    let successes: [Int] = probabilities.map { probability in
-      guard let index = cumulativeProbs.firstIndex(where: { cumulativeProb in
-          return probability <= cumulativeProb
-        })
-      else {
-        fatalError("The produced probability was not in range [0, 1]. Please file a bug report.")
-      }
-
-      return index
-    }
-
-    return successes
+    var rng = Xoroshiro256StarStar()
+    let weights = (0 ... trials).map { pmf(x: $0) }
+    let table = WalkersAliasTable(weights: weights)
+    return (0 ..< numberOfElements).map { _ in table.sample(using: &rng) }
   }
 }

Sources/StatKit/WalkersAliasTable.swift

@@ -0,0 +1,55 @@
+internal struct WalkersAliasTable {
+    let prob: [Double]
+    let alias: [Int]
+
+    init(weights: [Double]) {
+        let n = weights.count
+        let sum = weights.reduce(0, +)
+        // Normalise so probabilities sum to n (one per column)
+        var scaled = weights.map { $0 * Double(n) / sum }
+
+        var prob  = [Double](repeating: 0, count: n)
+        var alias = [Int](repeating: 0, count: n)
+
+        var small = [Int]()
+        var large = [Int]()
+
+        for i in 0 ..< n {
+            if scaled[i] < 1.0 {
+                small.append(i)
+            } else {
+                large.append(i)
+            }
+        }
+
+        while !small.isEmpty && !large.isEmpty {
+            let s = small.removeLast()
+            let l = large.removeLast()
+
+            prob[s]  = scaled[s]
+            alias[s] = l
+
+            // Donate the leftover from l to fill s's column
+            scaled[l] -= (1.0 - scaled[s])
+
+            if scaled[l] < 1.0 {
+                small.append(l)
+            } else {
+                large.append(l)
+            }
+        }
+
+        // Numerical rounding — anything left is effectively probability 1
+        for i in small + large { prob[i] = 1.0 }
+
+        self.prob  = prob
+        self.alias = alias
+    }
+
+    func sample(using rng: inout some RandomNumberGenerator) -> Int {
+        let n = prob.count
+        let i = Int.random(in: 0 ..< n, using: &rng)
+        let coin = Double.random(in: 0 ..< 1, using: &rng)
+        return coin < prob[i] ? i : alias[i]
+    }
+}

Tests/StatKitTests/Descriptive Statistics Tests/Distribution Tests/Discrete/BinomialDistributionTests.swift

@@ -105,8 +105,8 @@ struct BinomialDistributionTests {
 
   @Test("Sampling from a distribution returns correct proportions")
   func testSampling() async throws {
-    let numberOfSamples = 4000000
-    let numberOfTrials = 5
+    let numberOfSamples = 1000000
+    let numberOfTrials = 50
     let distribution = BinomialDistribution(probability: 0.7, trials: numberOfTrials)
     let samples = distribution.sample(numberOfSamples)
 
@@ -118,7 +118,7 @@ struct BinomialDistributionTests {
     let testRange =  0 ... numberOfTrials
 
     for successes in testRange {
-      #expect(proportions[successes, default: -1.0].isApproximatelyEqual(to: distribution.pmf(x: successes), absoluteTolerance: 1e-3))
+      #expect(proportions[successes, default: 0.0].isApproximatelyEqual(to: distribution.pmf(x: successes), absoluteTolerance: 1e-3))
     }
   }
 }