Replaced randga() with rgamma() (closes #9)

This commit is contained in:
Waldir Leoncio 2025-01-15 11:07:44 +01:00
parent 371ec35a06
commit 57dcff2209
5 changed files with 3 additions and 91 deletions

View file

@ -36,6 +36,7 @@ importFrom(matlab2r,times)
importFrom(matlab2r,uiputfile) importFrom(matlab2r,uiputfile)
importFrom(matlab2r,zeros) importFrom(matlab2r,zeros)
importFrom(methods,is) importFrom(methods,is)
importFrom(stats,rgamma)
importFrom(stats,runif) importFrom(stats,runif)
importFrom(stats,sd) importFrom(stats,sd)
importFrom(utils,read.delim) importFrom(utils,read.delim)

View file

@ -3,11 +3,11 @@
#' @examples rBAPS:::randdir(matrix(c(10, 30, 60), 3), 3) #' @examples rBAPS:::randdir(matrix(c(10, 30, 60), 3), 3)
#' @param counts shape parameter #' @param counts shape parameter
#' @param nc number of rows on output #' @param nc number of rows on output
#' @seealso randga #' @importFrom stats rgamma
randdir <- function(counts, nc) { randdir <- function(counts, nc) {
svar <- zeros(nc, 1) svar <- zeros(nc, 1)
for (i in 1:nc) { for (i in 1:nc) {
svar[i, 1] <- randga(counts[i, 1], 1) svar[i, 1] <- rgamma(1L, shape = counts[i, 1], rate = 1)
} }
svar <- svar / sum(svar) svar <- svar / sum(svar)
return(svar) return(svar)

View file

@ -1,61 +0,0 @@
#' @title Generates random number from a Gamma distribution
#' @description Generates one random number from shape parameter a and rate
#' parameter b
#' @param a shape
#' @param b rate
#' @return One realization of Gamma(a, b)
#' @details The generated random variable has mean a / b. It will be
#' positively-skewed for small values, but converges to a symmetric distribution
#' for very large numbers of a and b.
randga <- function(a, b) {
flag <- 0
if (a > 1) {
c1 <- a - 1
c2 <- (a - (1 / (6 * a))) / c1
c3 <- 2 / c1
c4 <- c3 + 2
c5 <- 1 / sqrt(a)
U1 <- 1
while (flag == 0) {
if (a <= 2.5) {
U1 <- rand()
U2 <- rand()
} else {
while (!(U1 > 0 & U1 < 1)) {
U1 <- rand()
U2 <- rand()
U1 <- U2 + c5 * (1 - 1.86 * U1)
}
}
W <- c2 * U2 / U1
if (c3 * U1 + W + (1 / W) <= c4) {
flag <- 1
g <- c1 * W / b
} else if (c3 * log(U1) - log(W) + W < 1) {
flag <- 1
g <- c1 * W / b
} else {
U1 <- -1
}
}
} else if (a == 1) {
g <- sum(-(1 / b) * log(rand(a, 1)))
} else {
while (flag == 0) {
U <- rand(2, 1)
if (U[1] > exp(1) / (a + exp(1))) {
g <- -log(((a + exp(1)) * (1 - U[1])) / (a * exp(1)))
if (U[2] <= g^(a - 1)) {
flag <- 1
}
} else {
g <- ((a + exp(1)) * U[1] / ((exp(1))^(1 / a)))
if (U[2] <= exp(-g)) {
flag <- 1
}
}
}
g <- g / b
}
return(g)
}

View file

@ -20,6 +20,3 @@ Generates random numbers
\examples{ \examples{
rBAPS:::randdir(matrix(c(10, 30, 60), 3), 3) rBAPS:::randdir(matrix(c(10, 30, 60), 3), 3)
} }
\seealso{
randga
}

View file

@ -1,25 +0,0 @@
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/randga.R
\name{randga}
\alias{randga}
\title{Generates random number from a Gamma distribution}
\usage{
randga(a, b)
}
\arguments{
\item{a}{shape}
\item{b}{rate}
}
\value{
One realization of Gamma(a, b)
}
\description{
Generates one random number from shape parameter a and rate
parameter b
}
\details{
The generated random variable has mean a / b. It will be
positively-skewed for small values, but converges to a symmetric distribution
for very large numbers of a and b.
}