diff --git a/NAMESPACE b/NAMESPACE
index 4b0f262..d49a2b7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,4 +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/rand.R b/R/rand.R
new file mode 100644
index 0000000..9e23039
--- /dev/null
+++ b/R/rand.R
@@ -0,0 +1,9 @@
+#' @title Generate matrix with U(0, 1) trials
+#' @description Imitates the behavior of `rand()` on Matlab
+#' @param r number of rows of output matrix
+#' @param c number of columns of output matrix
+#' @return \eqn{r \times c} matrix with random trials from a standard uniform distribution.
+#' @importFrom stats runif
+rand <- function(r = 1, c = 1) {
+    matrix(runif(r * c), r, c)
+}
\ No newline at end of file
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/TODO.md b/TODO.md
index f8522f6..d8c50d2 100644
--- a/TODO.md
+++ b/TODO.md
@@ -18,6 +18,8 @@ The list below contains non-essential but nice-to-have tasks for the next stable
 
 - [ ] Uniformize capitalization of function names
 - [ ] Uniformize capitalization of function arguments
+- [ ] Implement optimizations from `fastbaps`
+- [ ] Replace redundant functions (ex.: `randga`)
 
 # Known pitfalls
 
diff --git a/man/rand.Rd b/man/rand.Rd
new file mode 100644
index 0000000..131a43c
--- /dev/null
+++ b/man/rand.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rand.R
+\name{rand}
+\alias{rand}
+\title{Generate matrix with U(0, 1) trials}
+\usage{
+rand(r = 1, c = 1)
+}
+\arguments{
+\item{r}{number of rows of output matrix}
+
+\item{c}{number of columns of output matrix}
+}
+\value{
+\eqn{r \times c} matrix with random trials from a standard uniform distribution.
+}
+\description{
+Imitates the behavior of `rand()` on Matlab
+}
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.
+}