From d77ec891c50990d0dd38f7a6be448fda045e3359 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <waldir.leoncio@gmail.com>
Date: Tue, 17 Dec 2019 15:34:59 +0100
Subject: [PATCH] Added randga function

---
 NAMESPACE     |  1 +
 R/admix1.R    | 48 -----------------------------------------
 R/randga.R    | 59 +++++++++++++++++++++++++++++++++++++++++++++++++++
 man/randga.Rd | 22 +++++++++++++++++++
 4 files changed, 82 insertions(+), 48 deletions(-)
 create mode 100644 R/randga.R
 create mode 100644 man/randga.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 55f7c5f..d49a2b7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,5 +5,6 @@ export(calculatePopLogml)
 export(computeRows)
 export(learn_simple_partition)
 export(ownNum2Str)
+export(randga)
 export(repmat)
 importFrom(stats,runif)
diff --git a/R/admix1.R b/R/admix1.R
index 4453b16..e302f4f 100644
--- a/R/admix1.R
+++ b/R/admix1.R
@@ -695,54 +695,6 @@ admix1 <- function(tietue) {
 #     end;        
 # end;
 
-# %-------------------------------------------------------
-
-# function g=randga(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); 
-# end %while 
-# end %if
-# W = c2*U2/U1; 
-# if c3*U1+W+(1/W)<=c4, 
-# flag = 1; 
-# g = c1*W/b; 
-# elseif c3*log(U1)-log(W)+W<1, 
-# flag = 1; 
-# g = c1*W/b; 
-# else 
-# U1=-1; 
-# end %if 
-# end %while flag
-# elseif 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; 
-# end %if 
-# else 
-# g = ((a+exp(1))*U(1)/((exp(1))^(1/a))); 
-# if U(2)<=exp(-g), 
-# flag = 1; 
-# end %if 
-# end %if 
-# end %while flag 
-# g=g/b; 
-# end %if;
-
-
 # %-------------------------------------------------
 
 # function svar=randdir(counts,nc)
diff --git a/R/randga.R b/R/randga.R
new file mode 100644
index 0000000..2c721ad
--- /dev/null
+++ b/R/randga.R
@@ -0,0 +1,59 @@
+#' @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.
+#' @export
+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)
+}
\ No newline at end of file
diff --git a/man/randga.Rd b/man/randga.Rd
new file mode 100644
index 0000000..704a7a6
--- /dev/null
+++ b/man/randga.Rd
@@ -0,0 +1,22 @@
+% 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.
+}