Add Gamma distribution

.github/workflows/CodeStyle.yml

@@ -13,6 +13,6 @@ jobs:
     steps:
       - uses: actions/checkout@v1
       - name: Swiftlint
-        uses: norio-nomura/action-swiftlint@3.1.0
+        uses: norio-nomura/action-swiftlint@3.2.1
         with:
           args: --strict

.github/workflows/TestSuite.yml

@@ -5,22 +5,26 @@ on:
     branches: [ master ]
   push:
     branches: [ master ]
+
+env:
+  DEVELOPER_DIR: /Applications/Xcode_14.0.app
 
 jobs:
   macOS:
-    runs-on: macos-latest
+    runs-on: macos-12
 
     steps:
-    - name: Checkout
-      uses: actions/checkout@v2
+    - uses: actions/checkout@v3
     - name: Run macOS Tests
       run: swift test -c release --parallel
 
   Linux:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-20.04
 
     steps:
-      - name: Checkout
-        uses: actions/checkout@v2
-      - name: Run Linux Tests
-        run: swift test -c release --parallel
+    - uses: actions/checkout@v3
+    - uses: swift-actions/setup-swift@v1
+      with:
+        swift-version: '5.7'
+    - name: Run Linux Tests
+      run: swift test -c release --parallel

Package.swift

@@ -1,4 +1,4 @@
-// swift-tools-version:5.5
+// swift-tools-version:5.7
 
 import PackageDescription
 

README.md

@@ -19,7 +19,7 @@
 </p>
 
 <p align="center">
-    <img src="https://img.shields.io/badge/swift-5.5-blueviolet.svg" />
+    <img src="https://img.shields.io/badge/swift-5.7-blueviolet.svg" />
     <img src="https://img.shields.io/badge/swift pm-compatible-blueviolet.svg?style=flat" alt="Swift PM Compatible" />
 </p>
 

Sources/StatKit/ContinuedFraction.swift

@@ -0,0 +1,60 @@
+import RealModule
+
+/// A protocol outlining the functionality that a continued fraction generator must provide.
+internal protocol ContinuedFractionGenerator {
+  associatedtype NumericType: Real
+  /// A function generating parameters `alpha_n` for iterations in [1, +.infinity].
+  func alpha(for iteration: Int) -> NumericType
+  /// A function generating parameters `beta_n` for iterations in [1, +.infinity].
+  func beta(for iteration: Int) -> NumericType
+}
+
+/// Computes an estimate of a continued fraction, using the modified Lentz' algorithm.
+/// - parameter generator: A generator for the sequences of alpha's and beta's used in the computation.
+/// - returns: An estimate of the provided continued fraction.
+internal func continuedFraction<RealType, Generator>(generator: Generator) -> RealType
+where RealType: BinaryFloatingPoint,
+      Generator: ContinuedFractionGenerator,
+      Generator.NumericType == RealType
+{
+
+  let machineError = RealType.ulpOfOne
+  let ulpMultiplier = RealType(10)
+  let tiny = RealType.leastNonzeroMagnitude + .leastNonzeroMagnitude.ulp * ulpMultiplier
+
+  let a1 = generator.alpha(for: 1)
+  let b1 = generator.beta(for: 1)
+  let one = RealType(1)
+
+  var f = a1 / b1
+  var C = a1 / tiny
+  var D = one / b1
+  var iteration = 2
+  var delta = RealType.zero
+
+  repeat {
+    let an = generator.alpha(for: iteration)
+    let bn = generator.beta(for: iteration)
+
+    D = D * an + bn
+
+    if D == 0 {
+      D = tiny
+    }
+
+    C = bn + an / C
+
+    if C == 0 {
+      C = tiny
+    }
+
+    D = 1 / D
+    delta = C * D
+    f *= delta
+
+    iteration += 1
+
+  } while abs(delta - one) >= machineError
+
+  return f
+}

Sources/StatKit/Descriptive Statistics/Distribution/Continuous/GammaDistribution.swift

@@ -0,0 +1,187 @@
+/// A type modelling a Gamma Distribution.
+public struct GammaDistribution: ContinuousDistribution, UnivariateDistribution {
+  /// The Gamma distribution shape parameter.
+  public let alpha: Double
+
+  /// The Gamma distribution rate parameter.
+  public let beta: Double
+
+  /// Creates a Gamma Distribution with a specified alpha and beta parameters.
+  /// - parameter alpha: The distribution shape parameter.
+  /// - parameter beta: The distribution rate parameter.
+  public init(alpha: Double, beta: Double) {
+    precondition(0 < alpha, "The alpha parameter needs to be greater than 0.")
+    precondition(0 < beta, "The beta parameter needs to be greater than 0.")
+
+    self.alpha = alpha
+    self.beta = beta
+  }
+
+  /// Creates a Gamma Distribution with a specified shape and rate parameters.
+  /// - parameter shape: The distribution shape parameter.
+  /// - parameter rate: The distribution rate parameter.
+  public init(shape: Double, rate: Double) {
+    precondition(0 < shape, "The shape parameter needs to be greater than 0.")
+    precondition(0 < rate, "The rate parameter needs to be greater than 0.")
+
+    self.alpha = shape
+    self.beta = rate
+  }
+
+  public func pdf(x: Double, logarithmic: Bool = false) -> Double {
+    guard 0 < x else {
+      return logarithmic ? -.infinity : .zero
+    }
+
+    let logNumerator = (alpha - 1) * .log(x) - x / beta
+    let logDenominator = Double.logGamma(alpha) + alpha * .log(beta)
+    let logProbability = logNumerator - logDenominator
+
+    return logarithmic ? logProbability : .exp(logProbability)
+  }
+
+  public var mean: Double {
+    return alpha * beta
+  }
+
+  public var variance: Double {
+    return alpha * .pow(beta, 2)
+  }
+
+  public var skewness: Double {
+    return 2 / alpha.squareRoot()
+  }
+
+  public var kurtosis: Double {
+    return 6 / alpha + 3
+  }
+
+  public func cdf(x: Double, logarithmic: Bool = false) -> Double {
+    guard 0 < x else {
+      return logarithmic ? -.infinity : .zero
+    }
+    return regularizedIncompleteGamma(lower: .zero, upper: beta * x, alpha: alpha, logarithmic: logarithmic)
+  }
+
+  public func sample() -> Double {
+    var uniformGenerator = Xoroshiro256StarStar()
+
+    switch alpha {
+    case 0..<1:
+      return kunduGuptaSampling(using: &uniformGenerator)
+
+    default:
+      return martinoLuengoSampling(using: &uniformGenerator)
+    }
+  }
+
+  public func sample(_ numberOfElements: Int) -> [Double] {
+    precondition(0 < numberOfElements, "The requested number of samples need to be greater than 0.")
+    var uniformGenerator = Xoroshiro256StarStar()
+    let sampleRange = 1 ... numberOfElements
+
+    switch alpha {
+    case 0..<1:
+      return sampleRange.map { _ in kunduGuptaSampling(using: &uniformGenerator) }
+
+    default:
+      return sampleRange.map { _ in martinoLuengoSampling(using: &uniformGenerator) }
+    }
+  }
+
+  /// Computes a sample of a Gamma distribution with `0 < alpha < 1`.
+  /// - parameter generator: A uniform random number generator to be used in the computation.
+  /// - returns: A Gamma distribution sample.
+  ///
+  /// `alpha` must be in range (0, 1).
+  ///
+  /// Inspired by and based on work by Kundu and Gupta in their paper
+  /// '**A Convenient Way of Generating Gamma Random Variables Using Generalized Exponential Distribution**'.
+  private func kunduGuptaSampling(using generator: inout some RandomNumberGenerator) -> Double {
+    let d = 1.0334 - 0.0766 * .exp(2.2942 * alpha)
+
+    let firstFactorA = Double.pow(2.0, alpha)
+    let secondFactorA = Double.pow(1 - .exp(-d / 2), alpha)
+    let a = firstFactorA * secondFactorA
+
+    let b = alpha * .pow(d, alpha - 1) * .exp(-d)
+
+    let c = a + b
+
+    for _ in 1 ... 50 {
+      let U = Double.random(in: Double.leastNonzeroMagnitude ..< 1, using: &generator)
+
+      let X: Double
+
+      if U <= a / c {
+        X = -2 * .log(1 - .pow(c * U, 1 / alpha) / 2)
+      } else {
+        let numerator = c * (1 - U)
+        let denominator = alpha * .pow(d, alpha - 1)
+        X = -.log(numerator / denominator)
+      }
+
+      let V = Double.random(in: Double.leastNonzeroMagnitude ..< 1, using: &generator)
+      let threshold: Double
+
+      if X <= d {
+        let thresholdNumerator = .pow(X, alpha - 1) * .exp(-X / 2)
+        let thresholdDenominator = .pow(2, alpha - 1) * .pow(1 - .exp(-X / 2), alpha - 1)
+        threshold = thresholdNumerator / thresholdDenominator
+      } else {
+        threshold = .pow(d / X, 1 - alpha)
+      }
+
+      if V <= threshold {
+        return X
+      }
+    }
+
+    preconditionFailure("Failed to generate Gamma variate within a reasonable amount of time.")
+  }
+
+  /// Computes a sample of a Gamma distribution with `alpha >= 1`.
+  /// - parameter generator: A uniform random number generator to be used in the computation.
+  /// - returns: A Gamma distribution sample.
+  ///
+  /// `alpha` must be in range [1, +infinity].
+  /// This method uses rejection sampling to generate proposal samples from the Gamma distribution.
+  ///
+  /// Inspired by and based on work by Martino and Luengo in their paper
+  /// '**Extremely Efficient Generation of Gamma Random Variables for `alpha` >= 1**'.
+  private func martinoLuengoSampling(using generator: inout some RandomNumberGenerator) -> Double {
+    let proposalAlpha = alpha.rounded(.down)
+
+    let proposalBeta: Double
+    let proposalK: Double
+
+    switch alpha {
+    case 1..<2:
+      proposalBeta = beta / alpha
+      proposalK = .exp(1 - alpha) * .pow(alpha / beta, alpha - 1)
+
+    default:
+      proposalBeta = beta * (proposalAlpha - 1) / (alpha - 1)
+      proposalK = .exp(proposalAlpha - alpha) * .pow((alpha - 1) / beta, alpha - proposalAlpha)
+    }
+
+    for _ in 1 ... 50 {
+      let uniformProduct = (1 ... Int(proposalAlpha))
+        .map { _ in
+          Double.random(in: 0 ... 1, using: &generator)
+        }
+        .reduce(into: 1, *=)
+
+      let proposal = -.log(uniformProduct) / proposalBeta
+      let random = Double.random(in: 0 ... 1, using: &generator)
+      let p = .pow(proposal, alpha - 1) * .exp(-beta * proposal)
+      let pi = proposalK * .pow(proposal, proposalAlpha - 1) * .exp(-proposalBeta * proposal)
+
+      if random <= p / pi {
+        return proposal
+      }
+    }
+
+    preconditionFailure("Failed to generate Gamma variate within a reasonable amount of time.")
+  }
+}

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

@@ -33,7 +33,7 @@ public struct PoissonDistribution: DiscreteDistribution, UnivariateDistribution
   }
 
   public func pmf(x: Int, logarithmic: Bool = false) -> Double {
-    guard 0 <= x else { return 0 }
+    guard 0 <= x else { return logarithmic ? -.infinity : 0 }
 
     let logPMF = .log(.pow(rate, x.realValue)) - rate - .logGamma(x.realValue + 1)
     return logarithmic ? logPMF : .exp(logPMF)