From 15b11c59a2d8423b91ffcd5c4466d9c37af07ae9 Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmes@ufpr.br>
Date: Mon, 26 Oct 2015 18:03:21 -0300
Subject: [PATCH] Scratch of the rcg function (random count Gamma).
---
R/gammacount.R | 73 ++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 73 insertions(+)
create mode 100644 R/gammacount.R
diff --git a/R/gammacount.R b/R/gammacount.R
new file mode 100644
index 0000000..b446ea8
--- /dev/null
+++ b/R/gammacount.R
@@ -0,0 +1,73 @@
+## Generate Gamma-Count random variables.
+
+n=1L; shape=1; rate=1
+n=1L; shape=1:4; rate=11:14
+n=4L; shape=1:4; rate=11:12
+size <- c(n=n, sh=length(shape), rt=length(rate))
+
+rgc <- function(n=1L, shape=1, rate=1){
+ rGammaCount <- function(shape=shape, rate=rate){
+ ## To save the sum of Gamma ramdom variables.
+ s <- 0
+ ## To keep each Gamma random number.
+ x <- numeric(0)
+ repeat {
+ ## Generate a Gamma.
+ u <- rgamma(n=1L, shape=shape, rate=rate)
+ ## Update the sum.
+ s <- s+u
+ if (s>=1) break
+ ## Append the new value.
+ x <- append(x, values=u)
+ }
+ return(length(x))
+ }
+ sizes <- c(n, length(shape), length(rate))
+ if (max(sizes[-1])>1L){
+ message(paste(
+ "`n` was ignored because some",
+ "parameter has length > 1."))
+ shapeRate <- cbind(shape=shape, rate=rate)
+ y <- mapply(FUN=rGammaCount,
+ shape=shapeRate[, "shape"],
+ rate=shapeRate[, "rate"])
+ } else {
+ y <- replicate(n, rGammaCount(shape=shape, rate=rate))
+ }
+ return(y)
+}
+
+L <- list(
+ y1=rgc(1000, shape=1, rate=1),
+ y2=rgc(1000, shape=5, rate=5),
+ y3=rgc(1000, shape=0.2, rate=0.2))
+lapply(L, FUN=function(x) c(mean(x), var(x)))
+
+##----------------------------------------------------------------------
+
+rgc <- function(n=1L, shape=1, rate=1){
+ ## To save the sum of Gamma ramdom variables.
+ s <- 0
+ ## To keep each Gamma random number.
+ x <- numeric(0)
+ repeat {
+ ## Generate a Gamma.
+ u <- rgamma(n=1L, shape=shape, rate=rate)
+ ## Update the sum.
+ s <- s+u
+ if (s>=n) break
+ ## Append the new value.
+ x <- append(x, values=u)
+ }
+ ## Tabulate (count) occurrences in each interval.
+ y <- tabulate(ceiling(cumsum(x)))
+ ## This is a Gamma Count random variable.
+ return(y)
+}
+
+L <- list(
+ y1=rgc(500, shape=1, rate=1),
+ y2=rgc(500, shape=5, rate=5),
+ y3=rgc(500, shape=0.2, rate=0.2))
+
+lapply(L, FUN=function(x) c(mean(x), var(x)))
--
GitLab