Added randga function

NAMESPACE

@@ -5,5 +5,6 @@ export(calculatePopLogml)
 export(computeRows)
 export(learn_simple_partition)
 export(ownNum2Str)
+export(randga)
 export(repmat)
 importFrom(stats,runif)

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)

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)
+}

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.
+}