From ffab5b7d685781160ec8def7295dbb78d0ec8217 Mon Sep 17 00:00:00 2001
From: Yongfu Liao <liao961120@gmail.com>
Date: Wed, 17 Jul 2024 14:49:22 +0800
Subject: [PATCH] New Gamma-Poisson (negative binomial) distribution

---
 R/distributions.R | 21 ++++++++++++++++++++-
 1 file changed, 20 insertions(+), 1 deletion(-)

diff --git a/R/distributions.R b/R/distributions.R
index 3b01c33..dc5a043 100644
--- a/R/distributions.R
+++ b/R/distributions.R
@@ -14,7 +14,7 @@ rbern = function(prob, n=length(prob)) rbinom( n=n, size=1, prob=prob )
 #' The parameters are set to generate numbers from a half-normal distribution
 #' by default.
 #'
-#' @param n Numbers of sample.
+#' @param n Number of samples.
 #' @param m The mean of the untruncated normal distribution.
 #' @param s The standard deviation of the untruncated normal distribution.
 #' @param lower The lower bound for truncation.
@@ -55,3 +55,22 @@ rordlogit = function( phi, kappa ) {
   })
 }
 
+
+
+#' Random Number Generator from Gamma-Poisson Distribution
+#'
+#' @param n Number of samples.
+#' @param mu Poisson mean (rate)
+#' @param scale Parameter determining the dispersion of the outcomes.
+#'        This parameter is equivalent to Stan's `neg_binomial_2(mu, phi)`'s `phi`.
+#'
+#' @details
+#' See Chapter 12.1.2 of Statistical Rethinking 2th Edition, which details the
+#' Gamma-Poisson distribution for modeling over-dispersed count data.
+rgampois2 = function (n , mu , scale) {
+    # Matches stan neg_binomial_2(mu, phi=scale)
+    p_p = scale / (scale + mu)
+    p_n = scale
+    rnbinom(n, size = p_n, prob = p_p)
+}
+