Initial public release

.Rbuildignore

@@ -0,0 +1,8 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^README\.Rmd$
+^\.github$
+^data-raw$
+^_pkgdown\.yml$
+^docs$
+^pkgdown$

.aspell/defaults.R

@@ -0,0 +1,5 @@
+Rd_files <- vignettes <- R_files <- description <-
+  list(encoding = "UTF-8",
+       language = "en",
+       dictionaries = c("en_stats", "pg"))
+

.aspell/words.pws

@@ -0,0 +1,4 @@
+personal_ws-1.1 en 3 UTF-8
+Balamuta
+PG
+Polya-Gamma

.github/.gitignore

@@ -0,0 +1 @@
+*.html

.github/workflows/R-CMD-check.yaml

@@ -0,0 +1,29 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+
+name: R-CMD-check
+
+jobs:
+  R-CMD-check:
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
+    steps:
+      - uses: actions/checkout@v3
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::rcmdcheck
+          needs: check
+
+      - uses: r-lib/actions/check-r-package@v2

.gitignore

@@ -0,0 +1,8 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+src/*.o
+src/*.so
+src/*.dll
+docs

DESCRIPTION

@@ -0,0 +1,27 @@
+Package: pg
+Title: Polya-Gamma Distribution Sampler
+Version: 0.2.4
+Authors@R: 
+    c(
+      person(
+        given = "James",
+        family = "Balamuta",
+        role = c("aut", "cre", "cph"),
+        email = "balamut2@illinois.edu"
+      )
+    )
+Description: Provides access to a high performant random distribution
+    sampler for the Polya-Gamma Distribution using either C++ headers for 
+    'Rcpp' or 'RcppArmadillo' and 'R'.
+License: GPL (>= 3)
+LinkingTo: 
+    Rcpp,
+    RcppArmadillo
+Imports: 
+    Rcpp
+Encoding: UTF-8
+LazyData: true
+Roxygen: list(markdown = TRUE)
+RoxygenNote: 7.2.3
+Suggests: 
+    testthat (>= 2.1.0)

NAMESPACE

@@ -0,0 +1,13 @@
+# Generated by roxygen2: do not edit by hand
+
+export(pg_mean)
+export(pg_var)
+export(rpg_devroye)
+export(rpg_gamma)
+export(rpg_hybrid)
+export(rpg_normal)
+export(rpg_scalar)
+export(rpg_sp)
+export(rpg_vector)
+importFrom(Rcpp,sourceCpp)
+useDynLib(pg, .registration = TRUE)

NEWS.md

@@ -0,0 +1,3 @@
+# pg 0.2.4
+
+- Initial CRAN submission.

R/RcppExports.R

@@ -0,0 +1,111 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#' Sample from the Polya-Gamma distribution PG(h, z)
+#'
+#' Chooses the most efficient implemented method to sample from a Polya-Gamma
+#' distribution. Details on algorithm selection presented below.
+#'
+#' @param n     The number of samples to taken from a PG(h, z). Used only by
+#'              the vector sampler.
+#' @param h     `integer` values corresponding to the "shape" parameter.
+#' @param z     `numeric` values corresponding to the "scale" parameter.
+#' @param trunc Truncation cut-off. Only used by the gamma sampler.
+#'
+#' @return
+#' A single `numeric` value.
+#'
+#' @details
+#' The following sampling cases are enabled:
+#'
+#' - `h > 170`: Normal approximation method
+#' - `h > 13`: Saddlepoint approximation method
+#' - `h = 1` or `h = 2`: Devroye method
+#' - `h > 0`: Sum of Gammas method.
+#' - `h < 0`: Result is automatically set to zero.
+#'
+#' @export
+#' @rdname rpg-sampler
+#'
+#' @examples
+#' # Fixed parameter distribution simulation ----
+#'
+#' ## Parameters  ----
+#' h = 1; z = .5
+#'
+#' ## Sample only one value  ----
+#' single_value = rpg_scalar(h, z)
+#' single_value
+#'
+#' ## Attempt distribution recovery  ----
+#' vector_of_pg_samples = rpg_vector(1e6, h, z)
+#'
+#' head(vector_of_pg_samples)
+#' length(vector_of_pg_samples)
+#'
+#' ## Obtain the empirical results   ----
+#' empirical_mean = mean(vector_of_pg_samples)
+#' empirical_var = var(vector_of_pg_samples)
+#'
+#' ## Take the theoretical values ----
+#' theoretical_mean = pg_mean(h, z)
+#' theoretical_var = pg_var(h, z)
+#'
+#' ## Form a comparison table ----
+#'
+#' # empirically sampled vs. theoretical values
+#' rbind(c(empirical_mean, theoretical_mean),
+#'       c(empirical_var, theoretical_var))
+#'
+#' # Varying distribution parameters ----
+#'
+#' ## Generate varying parameters ----
+#' u_h = 20:100
+#' u_z = 0.5*u_h
+#'
+#' ## Sample from varying parameters ----
+#' x = rpg_hybrid(u_h, u_z)
+rpg_scalar <- function(h, z) {
+    .Call(`_pg_rpg_scalar`, h, z)
+}
+
+#' @export
+#' @rdname rpg-sampler
+rpg_vector <- function(n, h, z) {
+    .Call(`_pg_rpg_vector`, n, h, z)
+}
+
+#' @export
+#' @rdname rpg-sampler
+rpg_hybrid <- function(h, z) {
+    .Call(`_pg_rpg_hybrid`, h, z)
+}
+
+rpg_scalar_loop <- function(h, z) {
+    .Call(`_pg_rpg_scalar_loop`, h, z)
+}
+
+#' @export
+#' @rdname rpg-sampler
+rpg_gamma <- function(h, z, trunc = 1000L) {
+    .Call(`_pg_rpg_gamma`, h, z, trunc)
+}
+
+#' @export
+#' @rdname rpg-sampler
+rpg_devroye <- function(h, z) {
+    .Call(`_pg_rpg_devroye`, h, z)
+}
+
+#' @export
+#' @rdname rpg-sampler
+rpg_sp <- function(h, z) {
+    .Call(`_pg_rpg_sp`, h, z)
+}
+
+#' @export
+#' @rdname rpg-sampler
+rpg_normal <- function(h, z) {
+    .Call(`_pg_rpg_normal`, h, z)
+}
+

R/pg-package.R

@@ -0,0 +1,10 @@
+#' @keywords internal
+"_PACKAGE"
+
+# The following block is used by usethis to automatically manage
+# roxygen namespace tags. Modify with care!
+## usethis namespace: start
+#' @useDynLib pg, .registration = TRUE
+#' @importFrom Rcpp sourceCpp
+## usethis namespace: end
+NULL

R/pg_theoretical.R

@@ -0,0 +1,46 @@
+#' Theoretical Polya-Gamma Distribution's Mean and Variance
+#'
+#' Compute the theoretical mean and variance for a Polya-Gamma variable.
+#'
+#' @param h A single `integer` value corresponding to the "shape" parameter.
+#' @param z A single `numeric` value corresponding to the "scale" parameter.
+#'
+#' @return
+#' Either the theoretical mean or theoretical variance for a Polya-Gamma
+#' distribution.
+#'
+#' @export
+#' @rdname theoretical-pg
+#' @examples
+#' # Fixed parameter distribution simulation ----
+#'
+#' ## Parameters  ----
+#' h = 1; z = .5
+#' ## Attempt distribution recovery  ----
+#' vector_of_pg_samples = rpg_vector(1e6, h, z)
+#'
+#' head(vector_of_pg_samples)
+#' length(vector_of_pg_samples)
+#'
+#' ## Obtain the empirical results   ----
+#' empirical_mean = mean(vector_of_pg_samples)
+#' empirical_var = var(vector_of_pg_samples)
+#'
+#' ## Take the theoretical values ----
+#' theoretical_mean = pg_mean(h, z)
+#' theoretical_var = pg_var(h, z)
+#'
+#' ## Form a comparison table ----
+#'
+#' # empirically sampled vs. theoretical values
+#' rbind(c(empirical_mean, theoretical_mean),
+#'       c(empirical_var, theoretical_var))
+pg_mean = function(h, z) {
+  (0.5*h / z) * tanh(0.5 *z)
+}
+
+#' @export
+#' @rdname theoretical-pg
+pg_var = function(h, z) {
+  h / (4 * z ^ 3) * (sinh(z) - z) * (1 / cosh(0.5 * z) ^ 2)
+}