From 90af3d11a9dec8ff579a04f34ecbc6fb8fc0094e Mon Sep 17 00:00:00 2001 From: James J Balamuta Date: Mon, 17 Jul 2023 15:09:34 -0700 Subject: [PATCH] Initial public release --- .Rbuildignore | 8 + .aspell/defaults.R | 5 + .aspell/pg.rds | Bin 0 -> 72 bytes .aspell/words.pws | 4 + .github/.gitignore | 1 + .github/workflows/R-CMD-check.yaml | 29 +++ .gitignore | 8 + DESCRIPTION | 27 ++ NAMESPACE | 13 + NEWS.md | 3 + R/RcppExports.R | 111 ++++++++ R/pg-package.R | 10 + R/pg_theoretical.R | 46 ++++ README.Rmd | 167 ++++++++++++ README.md | 167 ++++++++++++ _pkgdown.yml | 4 + data-raw/custom-dictionary.R | 42 +++ inst/AUTHORS | 8 + inst/CITATION | 24 ++ inst/COPYRIGHTS | 17 ++ inst/WORDLIST | 3 + inst/include/pg.h | 63 +++++ inst/include/pg_inverse_gaussian_bones.h | 26 ++ inst/include/pg_inverse_gaussian_meat.h | 48 ++++ inst/include/pg_invert_y_bones.h | 80 ++++++ inst/include/pg_invert_y_meat.h | 122 +++++++++ inst/include/pg_polyagamma_bones.h | 69 +++++ inst/include/pg_polyagamma_meat.h | 266 +++++++++++++++++++ inst/include/pg_polyagamma_sp_bones.h | 66 +++++ inst/include/pg_polyagamma_sp_meat.h | 256 +++++++++++++++++++ inst/include/pg_r_rng_wrapper.h | 56 ++++ inst/include/pg_sampler_bones.h | 50 ++++ inst/include/pg_sampler_meat.h | 309 +++++++++++++++++++++++ inst/include/pg_truncated_gamma_bones.h | 33 +++ inst/include/pg_truncated_gamma_meat.h | 126 +++++++++ inst/include/pg_truncated_norm_bones.h | 50 ++++ inst/include/pg_truncated_norm_meat.h | 275 ++++++++++++++++++++ man/pg-package.Rd | 15 ++ man/rpg-sampler.Rd | 92 +++++++ man/theoretical-pg.Rd | 48 ++++ pg.Rproj | 21 ++ src/.gitignore | 3 + src/Makevars | 1 + src/Makevars.win | 1 + src/RcppExports.cpp | 129 ++++++++++ src/rpg.cpp | 125 +++++++++ tests/testthat.R | 4 + tests/testthat/test-saddlepoint.R | 49 ++++ tests/testthat/test-sample-control.R | 46 ++++ 49 files changed, 3126 insertions(+) create mode 100644 .Rbuildignore create mode 100644 .aspell/defaults.R create mode 100644 .aspell/pg.rds create mode 100644 .aspell/words.pws create mode 100644 .github/.gitignore create mode 100644 .github/workflows/R-CMD-check.yaml create mode 100644 .gitignore create mode 100644 DESCRIPTION create mode 100644 NAMESPACE create mode 100644 NEWS.md create mode 100644 R/RcppExports.R create mode 100644 R/pg-package.R create mode 100644 R/pg_theoretical.R create mode 100644 README.Rmd create mode 100644 README.md create mode 100644 _pkgdown.yml create mode 100644 data-raw/custom-dictionary.R create mode 100644 inst/AUTHORS create mode 100644 inst/CITATION create mode 100644 inst/COPYRIGHTS create mode 100644 inst/WORDLIST create mode 100644 inst/include/pg.h create mode 100644 inst/include/pg_inverse_gaussian_bones.h create mode 100644 inst/include/pg_inverse_gaussian_meat.h create mode 100644 inst/include/pg_invert_y_bones.h create mode 100644 inst/include/pg_invert_y_meat.h create mode 100644 inst/include/pg_polyagamma_bones.h create mode 100644 inst/include/pg_polyagamma_meat.h create mode 100644 inst/include/pg_polyagamma_sp_bones.h create mode 100644 inst/include/pg_polyagamma_sp_meat.h create mode 100644 inst/include/pg_r_rng_wrapper.h create mode 100644 inst/include/pg_sampler_bones.h create mode 100644 inst/include/pg_sampler_meat.h create mode 100644 inst/include/pg_truncated_gamma_bones.h create mode 100644 inst/include/pg_truncated_gamma_meat.h create mode 100644 inst/include/pg_truncated_norm_bones.h create mode 100644 inst/include/pg_truncated_norm_meat.h create mode 100644 man/pg-package.Rd create mode 100644 man/rpg-sampler.Rd create mode 100644 man/theoretical-pg.Rd create mode 100644 pg.Rproj create mode 100644 src/.gitignore create mode 100644 src/Makevars create mode 100644 src/Makevars.win create mode 100644 src/RcppExports.cpp create mode 100644 src/rpg.cpp create mode 100644 tests/testthat.R create mode 100644 tests/testthat/test-saddlepoint.R create mode 100644 tests/testthat/test-sample-control.R diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..40a07a8 --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,8 @@ +^.*\.Rproj$ +^\.Rproj\.user$ +^README\.Rmd$ +^\.github$ +^data-raw$ +^_pkgdown\.yml$ +^docs$ +^pkgdown$ diff --git a/.aspell/defaults.R b/.aspell/defaults.R new file mode 100644 index 0000000..064a80d --- /dev/null +++ b/.aspell/defaults.R @@ -0,0 +1,5 @@ +Rd_files <- vignettes <- R_files <- description <- + list(encoding = "UTF-8", + language = "en", + dictionaries = c("en_stats", "pg")) + diff --git a/.aspell/pg.rds b/.aspell/pg.rds new file mode 100644 index 0000000000000000000000000000000000000000..fadf6b01603ec1c8b8a83b26ef53018f64d0ee8c GIT binary patch literal 72 zcmb2|=3oE==I#ec2?+^F35jVb2}x{5k}M4~W;V7q3VRe=s6JD9Hf@@U(@!79cL8GM Y4t0L1i`G>q<>e_cXr3(m?F=*k07GXO+5i9m literal 0 HcmV?d00001 diff --git a/.aspell/words.pws b/.aspell/words.pws new file mode 100644 index 0000000..44100d9 --- /dev/null +++ b/.aspell/words.pws @@ -0,0 +1,4 @@ +personal_ws-1.1 en 3 UTF-8 +Balamuta +PG +Polya-Gamma diff --git a/.github/.gitignore b/.github/.gitignore new file mode 100644 index 0000000..2d19fc7 --- /dev/null +++ b/.github/.gitignore @@ -0,0 +1 @@ +*.html diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml new file mode 100644 index 0000000..f4b17a4 --- /dev/null +++ b/.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 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..45e67c6 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +.Rproj.user +.Rhistory +.RData +.Ruserdata +src/*.o +src/*.so +src/*.dll +docs diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..c810134 --- /dev/null +++ b/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) diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..30597be --- /dev/null +++ b/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) diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..e48189b --- /dev/null +++ b/NEWS.md @@ -0,0 +1,3 @@ +# pg 0.2.4 + +- Initial CRAN submission. diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..ce33e49 --- /dev/null +++ b/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) +} + diff --git a/R/pg-package.R b/R/pg-package.R new file mode 100644 index 0000000..4d9d176 --- /dev/null +++ b/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 diff --git a/R/pg_theoretical.R b/R/pg_theoretical.R new file mode 100644 index 0000000..9de66c7 --- /dev/null +++ b/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) +} diff --git a/README.Rmd b/README.Rmd new file mode 100644 index 0000000..ec8f3ee --- /dev/null +++ b/README.Rmd @@ -0,0 +1,167 @@ +--- +output: github_document +--- + + + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.path = "man/figures/README-", + out.width = "100%" +) +``` + +# pg + + +[![R-CMD-check](https://github.com/tmsalab/rcpppg/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/tmsalab/rcpppg/actions/workflows/R-CMD-check.yaml) +[![Package-License](http://img.shields.io/badge/license-GPL%20(%3E=3)-brightgreen.svg?style=flat)](http://www.gnu.org/licenses/gpl-3.0.html) + + +The goal of pg is to provide both _R_ and _C++_ header access to the Polya-Gamma +distribution sampling routine. + +## Installation + +You can install the development version of `pg` from +[GitHub](https://github.com/) with: + +```r +# install.packages("devtools") +devtools::install_github("tmsalab/pg") +``` + +## Usage + +Let `X` be a Polya-Gamma Distribution denoted by `PG(h, z)`, where `h` is the +"shape" parameter and `z` is the "scale" parameter. Presently, the following +sampling cases are enabled: + +- `h > 170`: Normal approximation method +- `h <= 170` and `h > 13`: Saddlepoint method +- `h = 1` or `h = 2`: Devroye method +- `h > 0`: Sum of gammas method. +- `h < 0`: Result is automatically set to zero. + +Not implemented: + +- `h <= 13` and `h > 1`: Alternative method (waiting for verification) + +The package structure allows for the sampling routines to be accessed +either via _C++_ or through _R_. The return type can be either a single +value or a vector. When repeat sampling is needed with the same `b` and `c`, +please use the vectorized sampler. + +### Sampling with C++ + +Using the sampling routine in _C++_ through a standalone `.cpp` file requires +either the `rpg_scalar_hybrid()`, `rpg_vector_hybrid()`, or `rpg_hybrid()` function to be +accessed in the `pg` _C++_ namespace. Each of these functions will automatically +select the appropriate algorithm based on criteria discussed previously. + +```cpp +#include +// [[Rcpp::depends(RcppArmadillo, pg)]] + +// [[Rcpp::export]] +double rpg_scalar(const double h, const double z) { + return pg::rpg_scalar_hybrid(h, z); +} + +// [[Rcpp::export]] +arma::vec rpg_hybrid(const arma::vec& h, const arma::vec& z) { + return pg::rpg_hybrid(h, z); +} + +// [[Rcpp::export]] +arma::vec rpg_vector(unsigned int n, const double h, const double z) { + return pg::rpg_vector_hybrid(n, h, z); +} +``` + +For use within an _R_ package, include a the `pg` package name in the +`DESCRIPTION` file. From there, include the `pg.h` header in a similar manner +to the stand-alone _C++_ example. + +``` +LinkingTo: + Rcpp, + RcppArmadillo + pg +``` + +### Sampling with R + +For use within an _R_ file, you can do: + +```{r} +# Number of observations to sample +n = 4 +# Select the PG(h, z) values +h = 1; z = 0.5 + +# Set a seed for reproducibility +set.seed(141) + +# Sample a single observation +pg::rpg_scalar(h, z) + +# Set a seed for reproducibility +set.seed(141) + +# Sample a vector of observations +pg::rpg_vector(n, h, z) +``` + + +## See also + +The following are useful resources regarding the Polya-Gamma distribution. + +- Papers + - "Bayesian Inference for Logistic Models Using Pólya–Gamma Latent Variables" + by Nicholas G. Polson, James G. Scott, and Jesse Windle (2013) + [doi:10.1080/01621459.2013.829001](https://www.tandfonline.com/doi/full/10.1080/01621459.2013.829001). Paper that invented the Polya-Gamma + - "Sampling Polya-Gamma random variates: alternate and approximate techniques" + by Jesse Windle, Nicholas G. Polson, and James G. Scott (2014) + . Provides an efficiency overview of the + different sampling approaches to sampling from a Polya-Gamma distribution. +- R Implementations: + - [`BayesLogit`](https://cran.r-project.org/package=BayesLogit) _R_ package by + Nicholas G. Polson, James G. Scott, and Jesse Windle. Provides the + main _C++_ class samplers contained used by the `pg` package. + - [`pgdraw`](https://cran.r-project.org/package=pgdraw) + by Daniel F. Schmidt and Enes Makalic. This package construction + relies heavily on free-floating functions and `Rcpp` data structures. + - [`bayesCL`](https://cran.r-project.org/package=bayesCL) + by Rok Cesnovar and Erik Strumbelj. This package boast a sampler that + is 100x faster through offloading of the computation onto a GPU. Though, + the package is not actively maintained. +- Support in other languages: + - `Python` has the [`pypolyagamma`](https://github.com/slinderman/pypolyagamma) + package by Scott Linderman. + - `Stan` [lacks an implementation](https://discourse.mc-stan.org/t/sampling-from-a-polya-gamma-distribution/8067) for the Polya-Gamma distribution since it relies on joint proposals + rather than full conditionals. + +## Author + +James Balamuta leaning heavily on work done in +[`BayesLogit`](https://cran.r-project.org/package=BayesLogit) _R_ package by +Nicholas G. Polson, James G. Scott, and Jesse Windle. + +## Citing the `pg` package + +To ensure future development of the package, please cite `pg` +package if used during an analysis or simulation study. Citation +information for the package may be acquired by using in *R*: + +``` r +citation("pg") +``` + +## License + +GPL (\>= 3) + diff --git a/README.md b/README.md new file mode 100644 index 0000000..d23d063 --- /dev/null +++ b/README.md @@ -0,0 +1,167 @@ + + + +# pg + + + +[![R-CMD-check](https://github.com/tmsalab/rcpppg/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/tmsalab/rcpppg/actions/workflows/R-CMD-check.yaml) +[![Package-License](http://img.shields.io/badge/license-GPL%20(%3E=3)-brightgreen.svg?style=flat)](http://www.gnu.org/licenses/gpl-3.0.html) + + +The goal of pg is to provide both *R* and *C++* header access to the +Polya-Gamma distribution sampling routine. + +## Installation + +You can install the development version of `pg` from +[GitHub](https://github.com/) with: + +``` r +# install.packages("devtools") +devtools::install_github("tmsalab/pg") +``` + +## Usage + +Let `X` be a Polya-Gamma Distribution denoted by `PG(h, z)`, where `h` +is the “shape” parameter and `z` is the “scale” parameter. Presently, +the following sampling cases are enabled: + +- `h > 170`: Normal approximation method +- `h <= 170` and `h > 13`: Saddlepoint method +- `h = 1` or `h = 2`: Devroye method +- `h > 0`: Sum of gammas method. +- `h < 0`: Result is automatically set to zero. + +Not implemented: + +- `h <= 13` and `h > 1`: Alternative method (waiting for verification) + +The package structure allows for the sampling routines to be accessed +either via *C++* or through *R*. The return type can be either a single +value or a vector. When repeat sampling is needed with the same `b` and +`c`, please use the vectorized sampler. + +### Sampling with C++ + +Using the sampling routine in *C++* through a standalone `.cpp` file +requires either the `rpg_scalar_hybrid()`, `rpg_vector_hybrid()`, or +`rpg_hybrid()` function to be accessed in the `pg` *C++* namespace. Each +of these functions will automatically select the appropriate algorithm +based on criteria discussed previously. + +``` cpp +#include +// [[Rcpp::depends(RcppArmadillo, pg)]] + +// [[Rcpp::export]] +double rpg_scalar(const double h, const double z) { + return pg::rpg_scalar_hybrid(h, z); +} + +// [[Rcpp::export]] +arma::vec rpg_hybrid(const arma::vec& h, const arma::vec& z) { + return pg::rpg_hybrid(h, z); +} + +// [[Rcpp::export]] +arma::vec rpg_vector(unsigned int n, const double h, const double z) { + return pg::rpg_vector_hybrid(n, h, z); +} +``` + +For use within an *R* package, include a the `pg` package name in the +`DESCRIPTION` file. From there, include the `pg.h` header in a similar +manner to the stand-alone *C++* example. + + LinkingTo: + Rcpp, + RcppArmadillo + pg + +### Sampling with R + +For use within an *R* file, you can do: + +``` r +# Number of observations to sample +n = 4 +# Select the PG(h, z) values +h = 1; z = 0.5 + +# Set a seed for reproducibility +set.seed(141) + +# Sample a single observation +pg::rpg_scalar(h, z) +#> [1] 0.05752942 + +# Set a seed for reproducibility +set.seed(141) + +# Sample a vector of observations +pg::rpg_vector(n, h, z) +#> [,1] +#> [1,] 0.05752942 +#> [2,] 0.38752679 +#> [3,] 0.38710433 +#> [4,] 0.18847913 +``` + +## See also + +The following are useful resources regarding the Polya-Gamma +distribution. + +- Papers + - “Bayesian Inference for Logistic Models Using Pólya–Gamma Latent + Variables” by Nicholas G. Polson, James G. Scott, and Jesse + Windle (2013) + [doi:10.1080/01621459.2013.829001](https://www.tandfonline.com/doi/full/10.1080/01621459.2013.829001). + Paper that invented the Polya-Gamma + - “Sampling Polya-Gamma random variates: alternate and approximate + techniques” by Jesse Windle, Nicholas G. Polson, and James G. + Scott (2014) . Provides an + efficiency overview of the different sampling approaches to sampling + from a Polya-Gamma distribution. +- R Implementations: + - [`BayesLogit`](https://cran.r-project.org/package=BayesLogit) *R* + package by Nicholas G. Polson, James G. Scott, and Jesse Windle. + Provides the main *C++* class samplers contained used by the `pg` + package. + - [`pgdraw`](https://cran.r-project.org/package=pgdraw) by Daniel F. + Schmidt and Enes Makalic. This package construction relies heavily + on free-floating functions and `Rcpp` data structures. + - [`bayesCL`](https://cran.r-project.org/package=bayesCL) by Rok + Cesnovar and Erik Strumbelj. This package boast a sampler that is + 100x faster through offloading of the computation onto a GPU. + Though, the package is not actively maintained. +- Support in other languages: + - `Python` has the + [`pypolyagamma`](https://github.com/slinderman/pypolyagamma) package + by Scott Linderman. + - `Stan` [lacks an + implementation](https://discourse.mc-stan.org/t/sampling-from-a-polya-gamma-distribution/8067) + for the Polya-Gamma distribution since it relies on joint proposals + rather than full conditionals. + +## Author + +James Balamuta leaning heavily on work done in +[`BayesLogit`](https://cran.r-project.org/package=BayesLogit) *R* +package by Nicholas G. Polson, James G. Scott, and Jesse Windle. + +## Citing the `pg` package + +To ensure future development of the package, please cite `pg` package if +used during an analysis or simulation study. Citation information for +the package may be acquired by using in *R*: + +``` r +citation("pg") +``` + +## License + +GPL (\>= 3) diff --git a/_pkgdown.yml b/_pkgdown.yml new file mode 100644 index 0000000..d71acfb --- /dev/null +++ b/_pkgdown.yml @@ -0,0 +1,4 @@ +url: ~ +template: + bootstrap: 5 + diff --git a/data-raw/custom-dictionary.R b/data-raw/custom-dictionary.R new file mode 100644 index 0000000..a0e8ac0 --- /dev/null +++ b/data-raw/custom-dictionary.R @@ -0,0 +1,42 @@ +custom_user_spelling_dictionary = function(custom_words, pkg) { + + # Check if package has aspell dictionary + aspell_location = ".aspell" + if(!dir.exists(aspell_location)) dir.create(aspell_location) + + # For newer versions of R, we create an R-level dictionary using an RDS file. + # This is read using + saveRDS(custom_words, file = file.path(aspell_location, paste0(pkg, ".rds")), version = 2) + + defaults_file_text = paste0('Rd_files <- vignettes <- R_files <- description <- + list(encoding = "UTF-8", + language = "en", + dictionaries = c("en_stats", "',pkg,'")) + ') + writeLines(defaults_file_text, file.path(aspell_location, "defaults.R")) + + # For the {spelling} package, write out a simple WORDLIST. + writeLines(custom_words, file.path("inst", "WORDLIST")) + + # Provide backward compatibility support for base R + # Note, must be run from package root. + # For more details, see: ?"aspell-utils" + utils::aspell_write_personal_dictionary_file( + custom_words, + out = file.path(aspell_location, "words.pws") + ) + + invisible(custom_words) +} + +author = c( + "Balamuta" +) +other = c( + "PG", + "Polya-Gamma" +) + +custom_words = c(author, other) + +custom_user_spelling_dictionary(custom_words, "pg") diff --git a/inst/AUTHORS b/inst/AUTHORS new file mode 100644 index 0000000..08bc05b --- /dev/null +++ b/inst/AUTHORS @@ -0,0 +1,8 @@ +The Polya-Gamma Distribution samplers were extracted from BayesLogit R package +with the following author and copyright information: + +Copyright: 2012-2019 Nicholas Polson, James Scott, Jesse Windle + +The pg package provides a refactorization of the algorithms to be self-contained +in a header-only library with additional C++ accessors. The refactorization +was done by James Joseph Balamuta. diff --git a/inst/CITATION b/inst/CITATION new file mode 100644 index 0000000..1a58d48 --- /dev/null +++ b/inst/CITATION @@ -0,0 +1,24 @@ +citHeader("To cite `pg` in publications please use use:") + +citEntry( + entry = "PhdThesis", + title = "Bayesian estimation of restricted latent class models: Extending priors, link functions, and structural models", + author = "James Joseph Balamuta", + school = "University of Illinois Urbana-Champaign", + doi = "10.1080/00273171.2021.1985949", + year = "2021", + textVersion = paste("Balamuta, J. J. (2021). Bayesian estimation of restricted latent class models: Extending priors, link functions, and structural models. University of Illinois Urbana-Champaign.") +) + +citEntry( + entry = "Unpublished", + title = "Bayesian inference for logistic models using Polya-Gamma latent variables", + author = "Nicholas G. Polson, James G. Scott, and Jesse Windle", + year = 2013, + url = "http://arxiv.org/abs/1205.0310", + note = "Most recent version: Jul. 2013.", + textVersion = + paste("Polson, N. G., Scott J. G., and Windle, J. (2013)", + "Bayesian inference for logistic models using Polya-Gamma latent variables.", + "URL .") +) diff --git a/inst/COPYRIGHTS b/inst/COPYRIGHTS new file mode 100644 index 0000000..7c31e8a --- /dev/null +++ b/inst/COPYRIGHTS @@ -0,0 +1,17 @@ +Overall license: +================ + + All of the work within the package is released under the + GNU GPL (>= 3). + +Details: +======== + +Files: * +Copyright: 2019 James Joseph Balamuta +License: GPL (>= 3) + +Files: inst/include/pg_* +Copyright: 2012-2019 Nicholas Polson, James Scott, Jesse Windle +License: GPL (>= 3) + diff --git a/inst/WORDLIST b/inst/WORDLIST new file mode 100644 index 0000000..679eccc --- /dev/null +++ b/inst/WORDLIST @@ -0,0 +1,3 @@ +Balamuta +PG +Polya-Gamma diff --git a/inst/include/pg.h b/inst/include/pg.h new file mode 100644 index 0000000..8b286f6 --- /dev/null +++ b/inst/include/pg.h @@ -0,0 +1,63 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_PG_H +#define PG_PG_H + +#include + +namespace pg { + +#ifndef USE_R +#define USE_R +#endif + +// Customized macros to work with RNGs +#include "pg_r_rng_wrapper.h" + +// Main polya-gamma sampling class +#include "pg_polyagamma_bones.h" +#include "pg_polyagamma_meat.h" + +// Utilities for saddlepoint estimation +#include "pg_inverse_gaussian_bones.h" +#include "pg_inverse_gaussian_meat.h" + +#include "pg_truncated_norm_bones.h" +#include "pg_truncated_norm_meat.h" + +#include "pg_truncated_gamma_bones.h" +#include "pg_truncated_gamma_meat.h" + +#include "pg_invert_y_bones.h" +#include "pg_invert_y_meat.h" + +// Saddlepoint-based polya-gamma sampling class +#include "pg_polyagamma_sp_bones.h" +#include "pg_polyagamma_sp_meat.h" + + +// Hybrid sampling definitions for C++ +#include "pg_sampler_bones.h" +#include "pg_sampler_meat.h" + +} // end namespace + +#endif + diff --git a/inst/include/pg_inverse_gaussian_bones.h b/inst/include/pg_inverse_gaussian_bones.h new file mode 100644 index 0000000..cae5d44 --- /dev/null +++ b/inst/include/pg_inverse_gaussian_bones.h @@ -0,0 +1,26 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_INVERSE_GAUSSIAN_BONES +#define PG_INVERSE_GAUSSIAN_BONES + +inline double igauss(double mu, double lambda); +inline double p_igauss(double x, double mu, double lambda); + +#endif diff --git a/inst/include/pg_inverse_gaussian_meat.h b/inst/include/pg_inverse_gaussian_meat.h new file mode 100644 index 0000000..78309d8 --- /dev/null +++ b/inst/include/pg_inverse_gaussian_meat.h @@ -0,0 +1,48 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_INVERSE_GAUSSIAN_MEAT +#define PG_INVERSE_GAUSSIAN_MEAT + +inline double igauss(double mu, double lambda) +{ + // See R code for specifics. + double mu2 = mu * mu; + double Y = r_norm(0.0, 1.0); + Y *= Y; + double W = mu + 0.5 * mu2 * Y / lambda; + double X = W - sqrt(W*W - mu2); + if (r_unif() > mu / (mu + X)) + X = mu2 / X; + return X; +} + + +inline double p_igauss(double x, double mu, double lambda) +{ + // z = 1 / mean + double z = 1 / mu; + double b = sqrt(lambda / x) * (x * z - 1); + double a = sqrt(lambda / x) * (x * z + 1) * -1.0; + // double y = p_norm(b, false) + exp(2 * lambda * z) * p_norm(a, false); + double y = p_norm(b, false) + exp(2 * lambda * z + p_norm(a, true)); + return y; +} + +#endif diff --git a/inst/include/pg_invert_y_bones.h b/inst/include/pg_invert_y_bones.h new file mode 100644 index 0000000..378507f --- /dev/null +++ b/inst/include/pg_invert_y_bones.h @@ -0,0 +1,80 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_INVERT_Y_BONES_H +#define PG_INVERT_Y_BONES_H + +const double tol = 1e-8; +const double IYPI = 3.141592653589793238462643383279502884197; + +inline double y_eval(double v); +inline void ydy_eval(double v, double* yp, double* dyp); + +inline double f_eval(double v, void* params); +inline double df_eval(double v, void* params); +inline void fdf_eval(double v, void* params, double* fp, double* dfp); + +inline double v_eval(double y, double tol=1e-9, int max_iter=1000); + +//------------------------------------------------------------------------------ + +const int grid_size = 81; +const double ygrid[] = { + 0.0625,0.06698584,0.07179365,0.07694653,0.08246924, + 0.08838835,0.09473229,0.1015315,0.1088188,0.1166291, + 0.125,0.1339717,0.1435873,0.1538931,0.1649385, + 0.1767767,0.1894646,0.2030631,0.2176376,0.2332582, + 0.25,0.2679434,0.2871746,0.3077861,0.329877, + 0.3535534,0.3789291,0.4061262,0.4352753,0.4665165, + 0.5,0.5358867,0.5743492,0.6155722,0.659754, + 0.7071068,0.7578583,0.8122524,0.8705506,0.933033, + 1,1.071773,1.148698,1.231144,1.319508, + 1.414214,1.515717,1.624505,1.741101,1.866066, + 2,2.143547,2.297397,2.462289,2.639016, + 2.828427,3.031433,3.24901,3.482202,3.732132, + 4,4.287094,4.594793,4.924578,5.278032, + 5.656854,6.062866,6.498019,6.964405,7.464264, + 8,8.574188,9.189587,9.849155,10.55606, + 11.31371,12.12573,12.99604,13.92881,14.92853, + 16}; + +const double vgrid[] = { + -256,-222.8609,-194.0117,-168.897,-147.0334, + -128,-111.4305,-97.00586,-84.4485,-73.51668, + -63.99997,-55.71516,-48.50276,-42.22387,-36.75755, + -31.99844,-27.85472,-24.24634,-21.10349,-18.36524, + -15.97843,-13.89663,-12.07937,-10.49137,-9.101928, + -7.884369,-6.815582,-5.875571,-5.047078,-4.315237, + -3.667256,-3.092143,-2.580459,-2.124095,-1.716085, + -1.350442,-1.022007,-0.7263359,-0.4595871,-0.2184366, + 0,0.1982309,0.3784427,0.5425468,0.6922181, + 0.828928,0.953973,1.068498,1.173516,1.269928, + 1.358533,1.440046,1.515105,1.584282,1.64809, + 1.706991,1.761401,1.811697,1.858218,1.901274, + 1.941143,1.978081,2.012318,2.044068,2.073521, + 2.100856,2.126234,2.149802,2.171696,2.192042, + 2.210954,2.228537,2.244889,2.260099,2.274249, + 2.287418,2.299673,2.311082,2.321703,2.331593, + 2.340804}; + +//------------------------------------------------------------------------------ + +//////////////////////////////////////////////////////////////////////////////// + +#endif diff --git a/inst/include/pg_invert_y_meat.h b/inst/include/pg_invert_y_meat.h new file mode 100644 index 0000000..18488e2 --- /dev/null +++ b/inst/include/pg_invert_y_meat.h @@ -0,0 +1,122 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_INVERT_Y_MEAT_H +#define PG_INVERT_Y_MEAT_H + +//------------------------------------------------------------------------------ + +inline double y_eval(double v) +{ + double y = 0.0; + double r = sqrt(fabs(v)); + if (v > tol) + y = tan(r) / r; + else if (v < -1*tol) + y = tanh(r) / r; + else + y = 1 + (1/3) * v + (2/15) * v * v + (17/315) * v * v * v; + return y; +} + +inline void ydy_eval(double v, double* yp, double* dyp) +{ + // double r = sqrt(fabs(v)); + + double y = y_eval(v); + *yp = y; + + if (fabs(v) >= tol) + *dyp = 0.5 * (y*y + (1-y) / v); + else + *dyp = 0.5 * (y*y - 1/3 - (2/15) * v); + +} + +inline double f_eval(double v, void * params) +{ + double y = *((double*) params); + return y_eval(v) - y; +} + +inline void fdf_eval(double v, void* params, double* fp, double* dfp) +{ + double y = *((double*)params); + ydy_eval(v, fp, dfp); + *fp -= y; +} + +inline double df_eval(double v, void * params) +{ + double f, df; + ydy_eval(v, &f, &df); + return df; +} + +inline double v_eval(double y, double tol, int max_iter) +{ + double ylower = ygrid[0]; + double yupper = ygrid[grid_size-1]; + + if (y < ylower) { + return -1. / (y*y); + } else if (y > yupper) { + double v = atan(0.5 * y * IYPI); + return v*v; + } + else if (y==1) return 0.0; + + double id = (log(y) / log(2.0) + 4.0) / 0.1; + // printf("y, id, y[id], v[id]: %g, %g, %g, %g\n", y, id, ygrid[(int)id], vgrid[(int)id]); + + // C++ default is truncate decimal portion. + int idlow = (int)id; + int idhigh = (int)id + 1; + double vl = vgrid[idlow]; // lower bound + double vh = vgrid[idhigh]; // upper bound + + int iter = 0; + double diff = tol + 1.0; + double vnew = vl; + double vold = vl; + double f0, f1; + + while (diff > tol && iter < max_iter) { + iter++; + vold = vnew; + fdf_eval(vold, &y, &f0, &f1); + vnew = vold - f0 / f1; + vnew = vnew > vh ? vh : vnew; + vnew = vnew < vl ? vl : vnew; + diff = fabs(vnew - vold); + // printf("iter: %i, v: %g, diff: %g\n", iter, vnew, diff); + } + + if (iter >= max_iter) { + #ifndef USE_R + fprintf(stderr, "InvertY.cpp, v_eval: reached max_iter: %i\n", iter); + #else + Rprintf("InvertY.cpp, v_eval: reached max_iter: %i\n", iter); + #endif + } + + return vnew; +} + +#endif diff --git a/inst/include/pg_polyagamma_bones.h b/inst/include/pg_polyagamma_bones.h new file mode 100644 index 0000000..bd33107 --- /dev/null +++ b/inst/include/pg_polyagamma_bones.h @@ -0,0 +1,69 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +// See for implementation details. + +#ifndef PG_POLYAGAMMA_BONES_H +#define PG_POLYAGAMMA_BONES_H + +// The numerical accuracy of __PI will affect your distribution. +const double __PI = 3.141592653589793238462643383279502884197; +const double HALFPISQ = 0.5 * __PI * __PI; +const double FOURPISQ = 4 * __PI * __PI; +const double __TRUNC = 0.64; +const double __TRUNC_RECIP = 1.0 / __TRUNC; + +class PolyaGamma +{ + + // For sum of Gammas. + int T; + std::vector bvec; + + public: + + // Constructors. + PolyaGamma(int trunc = 200); + + // Draw. + // double draw(double n, double z, RNG& r); + double draw(int n, double z); + double draw_sum_of_gammas(double n, double z); + double draw_like_devroye(double z); + + //void draw(MF x, double a, double z, RNG& r); + //void draw(MF x, MF a, MF z, RNG& r); + + // Utility. + void set_trunc(int trunc); + + // Helper. + double a(int n, double x); + double pigauss(double x, double Z); + double mass_texpon(double Z); + double rtigauss(double Z); + + static double jj_m1(double b, double z); + static double jj_m2(double b, double z); + static double pg_m1(double b, double z); + static double pg_m2(double b, double z); + +}; + +#endif diff --git a/inst/include/pg_polyagamma_meat.h b/inst/include/pg_polyagamma_meat.h new file mode 100644 index 0000000..52735cc --- /dev/null +++ b/inst/include/pg_polyagamma_meat.h @@ -0,0 +1,266 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_POLYAGAMMA_MEAT_H +#define PG_POLYAGAMMA_MEAT_H + +//////////////////////////////////////////////////////////////////////////////// + // Constructors // +//////////////////////////////////////////////////////////////////////////////// + +inline PolyaGamma::PolyaGamma(int trunc) : T(trunc), bvec(T) +{ + set_trunc(T); +} // PolyaGamma + +//////////////////////////////////////////////////////////////////////////////// + // Utility // +//////////////////////////////////////////////////////////////////////////////// + +inline void PolyaGamma::set_trunc(int trunc) +{ + + if (trunc < 1) { + #ifndef NTHROW + throw std::invalid_argument("PolyaGamma(int trunc): trunc < 1."); + #else + #ifndef USE_R + fprintf(stderr, "Invalid parameter: PolyaGamma(int trunc): trunc < 1. Setting trunc=1.\n"); + #else + Rprintf("Invalid parameter: PolyaGamma(int trunc): trunc < 1. Setting trunc=1.\n"); + #endif + trunc = 1; + #endif + } + + T = trunc; + bvec.resize(T); + + for(int k=0; k < T; ++k){ + // + since we start indexing at 0. + double d = ((double) k + 0.5); + bvec[k] = FOURPISQ * d * d; + } +} // set_trunc + +inline double PolyaGamma::a(int n, double x) +{ + double K = (n + 0.5) * __PI; + double y = 0; + if (x > __TRUNC) { + y = K * exp( -0.5 * K*K * x ); + } + else if (x > 0) { + double expnt = -1.5 * (log(0.5 * __PI) + log(x)) + log(K) - 2.0 * (n+0.5)*(n+0.5) / x; + y = exp(expnt); + // y = std::pow(0.5 * __PI * x, -1.5) * K * exp( -2.0 * (n+0.5)*(n+0.5) / x); + // ^- unstable for small x? + } + return y; +} + +inline double PolyaGamma::pigauss(double x, double Z) +{ + double b = sqrt(1.0 / x) * (x * Z - 1); + double a = sqrt(1.0 / x) * (x * Z + 1) * -1.0; + double y = p_norm(b, false) + exp(2 * Z) * p_norm(a, false); + return y; +} + +inline double PolyaGamma::mass_texpon(double Z) +{ + double t = __TRUNC; + + double fz = 0.125 * __PI*__PI + 0.5 * Z*Z; + double b = sqrt(1.0 / t) * (t * Z - 1); + double a = sqrt(1.0 / t) * (t * Z + 1) * -1.0; + + double x0 = log(fz) + fz * t; + double xb = x0 - Z + p_norm(b, true); + double xa = x0 + Z + p_norm(a, true); + + double qdivp = 4 / __PI * ( exp(xb) + exp(xa) ); + + return 1.0 / (1.0 + qdivp); +} + +inline double PolyaGamma::rtigauss(double Z) +{ + Z = fabs(Z); + double t = __TRUNC; + double X = t + 1.0; + if (__TRUNC_RECIP > Z) { // mu > t + double alpha = 0.0; + while (r_unif() > alpha) { + // X = t + 1.0; + // while (X > t) + // X = 1.0 / gamma_rate(0.5, 0.5); + // Slightly faster to use truncated normal. + double E1 = expon_rate(1.0); double E2 = expon_rate(1.0); + while ( E1*E1 > 2 * E2 / t) { + E1 = expon_rate(1.0); E2 = expon_rate(1.0); + } + X = 1 + E1 * t; + X = t / (X * X); + alpha = exp(-0.5 * Z*Z * X); + } + } + else { + double mu = 1.0 / Z; + while (X > t) { + double Y = r_norm(0., 1.); Y *= Y; + double half_mu = 0.5 * mu; + double mu_Y = mu * Y; + X = mu + half_mu * mu_Y - half_mu * sqrt(4 * mu_Y + mu_Y * mu_Y); + if (r_unif() > mu / (mu + X)) + X = mu*mu / X; + } + } + return X; +} + +//////////////////////////////////////////////////////////////////////////////// + // Sample // +//////////////////////////////////////////////////////////////////////////////// + +// double PolyaGamma::draw(double n, double z, RNG& r) +// { +// return draw_sum_of_gammas(n, z, r); +// } + +inline double PolyaGamma::draw(int n, double z) +{ + if (n < 1) { + #ifndef NTHROW + throw std::invalid_argument("PolyaGamma::draw: n < 1."); + #else + #ifndef USE_R + fprintf(stderr, "PolyaGamma::draw: n < 1. Set n = 1.\n"); + #else + Rprintf("PolyaGamma::draw: n < 1. Set n = 1.\n"); + #endif + n = 1; + #endif + } + double sum = 0.0; + for (int i = 0; i < n; ++i) + sum += draw_like_devroye(z); + return sum; +} // draw + +inline double PolyaGamma::draw_sum_of_gammas(double n, double z) +{ + double x = 0; + double kappa = z * z; + for(int k=0; k < T; ++k) + x += gamma_scale(n, 1.0) / (bvec[k] + kappa); + return 2.0 * x; +} // draw_sum_of_gammas + +inline double PolyaGamma::draw_like_devroye(double Z) +{ + // Change the parameter. + Z = fabs(Z) * 0.5; + + // Now sample 0.25 * J^*(1, Z := Z/2). + double fz = 0.125 * __PI*__PI + 0.5 * Z*Z; + // ... Problems with large Z? Try using q_over_p. + // double p = 0.5 * __PI * exp(-1.0 * fz * __TRUNC) / fz; + // double q = 2 * exp(-1.0 * Z) * pigauss(__TRUNC, Z); + + double X = 0.0; + double S = 1.0; + double Y = 0.0; + // int iter = 0; If you want to keep track of iterations. + + while (true) { + + // if (r_unif() < p/(p+q)) + if ( r_unif() < mass_texpon(Z) ) + X = __TRUNC + expon_rate(1.) / fz; + else + X = rtigauss(Z); + + S = a(0, X); + Y = r_unif() * S; + int n = 0; + bool go = true; + + // Cap the number of iterations? + while (go) { + + // Break infinite loop. Put first so it always checks n==0. + #ifdef USE_R + if (n % 1000 == 0) R_CheckUserInterrupt(); + #endif + + ++n; + if (n%2==1) { + S = S - a(n, X); + if ( Y<=S ) return 0.25 * X; + } + else { + S = S + a(n, X); + if ( Y>S ) go = false; + } + + } + // Need Y <= S in event that Y = S, e.g. when X = 0. + + } +} // draw_like_devroye + +//////////////////////////////////////////////////////////////////////////////// + // Static Members // +//////////////////////////////////////////////////////////////////////////////// + +inline double PolyaGamma::jj_m1(double b, double z) +{ + z = fabs(z); + double m1 = 0.0; + if (z > 1e-12) + m1 = b * tanh(z) / z; + else + m1 = b * (1 - (1.0/3) * std::pow(z,2) + (2.0/15) * std::pow(z,4) - (17.0/315) * std::pow(z,6)); + return m1; +} + +inline double PolyaGamma::jj_m2(double b, double z) +{ + z = fabs(z); + double m2 = 0.0; + if (z > 1e-12) + m2 = (b+1) * b * std::pow(tanh(z)/z,2) + b * ((tanh(z)-z)/std::pow(z,3)); + else + m2 = (b+1) * b * std::pow(1 - (1.0/3) * std::pow(z,2) + (2.0/15) * std::pow(z,4) - (17.0/315) * std::pow(z,6), 2) + + b * ((-1.0/3) + (2.0/15) * std::pow(z,2) - (17.0/315) * std::pow(z,4)); + return m2; +} + +inline double PolyaGamma::pg_m1(double b, double z) +{ + return jj_m1(b, 0.5 * z) * 0.25; +} + +inline double PolyaGamma::pg_m2(double b, double z) +{ + return jj_m2(b, 0.5 * z) * 0.0625; +} + +#endif diff --git a/inst/include/pg_polyagamma_sp_bones.h b/inst/include/pg_polyagamma_sp_bones.h new file mode 100644 index 0000000..9b25640 --- /dev/null +++ b/inst/include/pg_polyagamma_sp_bones.h @@ -0,0 +1,66 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_POLYAGAMMASP_BONES_H +#define PG_POLYAGAMMASP_BONES_H + +// Function and derivative. +struct FD { + double val; + double der; +}; + +struct Line { + double slope; + double icept; +}; + + +// PolyaGamma approximation by SP. +class PolyaGammaApproxSP +{ + +public: + + int draw(double& d, double h, double z, int max_iter=200); + +protected: + + // Helper. + + double w_left (double trunc, double h, double z); + double w_right(double trunc, double h, double z); + + void delta_func(double x, double z , FD& delta); + double phi_func (double x, double mid, FD& phi); + + double tangent_to_eta(double x, double z, double mid, Line& tl); + + double sp_approx(double x, double n, double z); + + double cos_rt(double v); + + // YV yv; + + double rtigauss(double mu, double lambda, double trunc); + double y_func(double v); // y = tan(sqrt(v)) / sqrt(v); + +}; + +#endif diff --git a/inst/include/pg_polyagamma_sp_meat.h b/inst/include/pg_polyagamma_sp_meat.h new file mode 100644 index 0000000..01ce38d --- /dev/null +++ b/inst/include/pg_polyagamma_sp_meat.h @@ -0,0 +1,256 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_POLYAGAMMASP_MEAT_H +#define PG_POLYAGAMMASP_MEAT_H + +inline double PolyaGammaApproxSP::rtigauss(double mu, double lambda, double trunc) +{ + // mu = fabs(mu); + double X = trunc + 1.0; + if (trunc < mu) { // mu > t + double alpha = 0.0; + while (r_unif() > alpha) { + X = rtinvchi2(lambda, trunc); + alpha = exp(-0.5 * lambda / (mu*mu) * X); + } + // printf("rtigauss, part i: %g\n", X); + } + else { + while (X > trunc) { + X = igauss(mu, lambda); + } + // printf("rtigauss, part ii: %g\n", X); + } + return X; +} + +inline double PolyaGammaApproxSP::y_func(double v) +{ + double tol = 1e-6; + double y = 0.0; + double r = sqrt(fabs(v)); + if (v > tol) + y = tan(r) / r; + else if (v < -1*tol) + y = tanh(r) / r; + else + y = 1 + (1/3) * v + (2/15) * v * v + (17/315) * v * v * v; + return y; +} + +inline double PolyaGammaApproxSP::cos_rt(double v) +{ + double y = 0.0; + double r = sqrt(fabs(v)); + if (v >= 0) + y = cos(r); + else + y = cosh(r); + return y; +} + +inline void PolyaGammaApproxSP::delta_func(double x, double mid, FD& delta) +{ + if (x >= mid) { + delta.val = log(x) - log(mid); + delta.der = 1.0 / x; + } + else { + delta.val = 0.5 * (1 - 1.0 / x) - 0.5 * (1 - 1.0 / mid); + delta.der = 0.5 / (x*x); + } +} + +inline double PolyaGammaApproxSP::phi_func(double x, double z, FD& phi) +{ + // double v = yv.v_func(x); + double v = v_eval(x); + double u = 0.5 * v; + double t = u + 0.5 * z*z; + + phi.val = log(cosh(fabs(z))) - log(cos_rt(v)) - t * x; + phi.der = -1.0 * t; + + return v; +} + +inline double PolyaGammaApproxSP::tangent_to_eta(double x, double z, double mid, Line& tl) +{ + FD phi, delta, eta; + double v; + + v = phi_func(x, z, phi); + delta_func(x, mid, delta); + + eta.val = phi.val - delta.val; + eta.der = phi.der - delta.der; + + // printf("v=%g\nphi=%g, phi.d=%g\ndelta=%g, delta.d=%g\neta=%g, eta.d=%g\n", + // v, phi.val, phi.der, delta.val, delta.der, eta.val, eta.der); + + tl.slope = eta.der; + tl.icept = eta.val - eta.der * x; + + return v; +} + +inline double PolyaGammaApproxSP::sp_approx(double x, double n, double z) +{ + // double v = yv.v_func(x); + double v = v_eval(x); + double u = 0.5 * v; + double z2 = z * z; + double t = u + 0.5 * z2; + // double m = y_func(-1 * z2); + + double phi = log(cosh(z)) - log(cos_rt(v)) - t * x; + + double K2 = 0.0; + if (fabs(v) >= 1e-6) + K2 = x*x + (1-x) / v; + else + K2 = x*x - 1/3 - (2/15) * v; + + double log_spa = 0.5 * log(0.5 * n / __PI) - 0.5 * log(K2) + n * phi; + return exp(log_spa); +} + +inline int PolyaGammaApproxSP::draw(double& d, double n, double z, int maxiter) +{ + if (n < 1) { + #ifndef USE_R + fprintf(stderr, "PolyaGammaApproxSP::draw: n must be >= 1.\n"); + #else + Rprintf("PolyaGammaApproxSP::draw: n must be >= 1.\n"); + #endif + return -1.; + } + + z = 0.5 * fabs(z); + + double xl = y_func(-1*z*z); // Mode of phi - Left point. + double md = xl * 1.1; // Mid point. + double xr = xl * 1.2; // Right point. + + // printf("xl, md, xr: %g, %g, %g\n", xl, md, xr); + + // Inflation constants + // double vmd = yv.v_func(md); + double vmd = v_eval(md); + double K2md = 0.0; + + if (fabs(vmd) >= 1e-6) + K2md = md*md + (1-md) / vmd; + else + K2md = md*md - 1/3 - (2/15) * vmd; + + double m2 = md * md; + double al = m2*md / K2md; + double ar = m2 / K2md; + + // printf("vmd, K2md, al, ar: %g, %g %g %g\n", vmd, K2md, al, ar); + + // Tangent lines info. + Line ll, lr; + tangent_to_eta(xl, z, md, ll); + tangent_to_eta(xr, z, md, lr); + + double rl = -1. * ll.slope; + double rr = -1. * lr.slope; + double il = ll.icept; + double ir = lr.icept; + + // printf("rl, rr, il, ir: %g, %g, %g, %g\n", rl, rr, il, ir); + + // Constants + double lcn = 0.5 * log(0.5 * n / __PI); + double rt2rl = sqrt(2 * rl); + + // printf("sqrt(rl): %g\n", rt2rl); + + // Weights + double wl, wr, wt, pl; + + + // // to cross-reference R script + // double term1, term2, term3, term4, term5; + // term1 = exp(0.5 * log(al)); + // term2 = exp(- n * rt2rl + n * il + 0.5 * n * 1./md); + // term3 = p_igauss(md, 1./rt2rl, n); + // printf("l terms 1-3: %g, %g, %g\n", term1, term2, term3); + + wl = exp(0.5 * log(al) - n * rt2rl + n * il + 0.5 * n * 1./md) * + p_igauss(md, 1./rt2rl, n); + + // // to cross-reference R script + // term1 = exp(0.5 * log(ar)); + // term2 = exp(lcn); + // term3 = exp(- n * log(n * rr) + n * ir - n * log(md)); + // term4 = exp(ln_gamma(n)); + // term5 = (1.0 - p_gamma_rate(md, n, n*rr, false)); + // printf("r terms 1-5: %g, %g, %g, %g, %g\n", term1, term2, term3, term4, term5); + + wr = exp(0.5 * log(ar) + lcn + (- n * log(n * rr) + n * ir - n * log(md) + ln_gamma(n)) ) * + (1.0 - p_gamma_rate(md, n, n*rr, false)); + // yv.upperIncompleteGamma(md, n, n*rr); + + // printf("wl, wr, lcn: %g, %g, %g\n", wl, wr, lcn); + + wt = wl + wr; + pl = wl / wt; + + // Sample + bool go = true; + int iter = 0; + double X = 2.0; + double F = 0.0; + + while(go && iter < maxiter) { + // Put first so check on first pass. + #ifdef USE_R + if (iter % 1000 == 0) R_CheckUserInterrupt(); + #endif + + iter++; + + double phi_ev; + if (r_unif() < pl) { + X = rtigauss(1./rt2rl, n, md); + phi_ev = n * (il - rl * X) + 0.5 * n * ((1.-1./X) - (1.-1./md)); + F = exp(0.5 * log(al) + lcn - 1.5 * log(X) + phi_ev); + } + else { + X = ltgamma(n, n*rr, md); + phi_ev = n * (ir - rr * X) + n * (log(X) - log(md)); + F = exp(0.5 * log(ar) + lcn + phi_ev) / X; + } + + double spa = sp_approx(X, n, z); + + if (F * r_unif() < spa) go = false; + + } + + // return n * 0.25 * X; + d = n * 0.25 * X; + return iter; +} + +#endif diff --git a/inst/include/pg_r_rng_wrapper.h b/inst/include/pg_r_rng_wrapper.h new file mode 100644 index 0000000..e54e30a --- /dev/null +++ b/inst/include/pg_r_rng_wrapper.h @@ -0,0 +1,56 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + + +// MAKE SURE YOU CALL GetRNGSeed() and PutRNGSeed() WHEN USING THESE FUNCTIONS!!! + +#ifndef SIMPLE_RNG_WRAPPER +#define SIMPLE_RNG_WRAPPER + +#define RCHECK 1000 + +const double SQRT2PI = 2.50662827; + +inline void check_R_interupt(int count) +{ + #ifdef USE_R + if (count % RCHECK == 0) R_CheckUserInterrupt(); + #endif +} + +// CDF +#define p_norm(X, USE_LOG) R::pnorm((X), 0.0, 1.0, true, (USE_LOG)) +#define p_gamma_scale(X, SHAPE, SCALE, USE_LOG) R::pgamma((X), (SHAPE), (SCALE), 1, (USE_LOG)) +#define p_gamma_rate(X, SHAPE, RATE, USE_LOG) R::pgamma((X), (SHAPE), (1./(RATE)), 1, (USE_LOG)) + +// Random variates +#define expon_mean(MEAN) R::rexp((MEAN)) +#define expon_rate(RATE) R::rexp((1./(RATE))) +#define r_unif() R::runif(0.0, 1.0) +#define r_norm(MEAN, SD) R::rnorm((MEAN), (SD)) +#define gamma_scale(SHAPE, SCALE) R::rgamma((SHAPE), (SCALE)) +#define gamma_rate(SHAPE, RATE) R::rgamma((SHAPE), (1./(RATE))) + +// Scientific functions +#define ln_gamma(X) R::lgammafn(X) + +// Output + + +#endif diff --git a/inst/include/pg_sampler_bones.h b/inst/include/pg_sampler_bones.h new file mode 100644 index 0000000..67d474a --- /dev/null +++ b/inst/include/pg_sampler_bones.h @@ -0,0 +1,50 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_SAMPLER_BONES_H +#define PG_SAMPLER_BONES_H + +template +inline T rpg_gamma(const T& shape_, const T& scale_, int trunc = 1000); + +// Ensure argument parity with other samplers by embedding trunc with a default of 1000. +template +inline T rpg_gamma_trunc(const T& shape_, const T& scale_); + +template +inline T rpg_devroye(const T& shape_, const T& scale_); + +//template +//inline T rpg_alt(const T& shape_, const T& scale_); + +template +inline T rpg_sp(const T& shape_, const T& scale_); + +template +inline T rpg_normal(const T& shape_, const T& scale_); + +inline arma::vec rpg_vector_hybrid( + unsigned int n, double shape, double scale); + +inline double rpg_scalar_hybrid(double shape, double scale); + +template +inline T rpg_hybrid(const T& shape_, const T& scale_); + +#endif diff --git a/inst/include/pg_sampler_meat.h b/inst/include/pg_sampler_meat.h new file mode 100644 index 0000000..ba7f045 --- /dev/null +++ b/inst/include/pg_sampler_meat.h @@ -0,0 +1,309 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_SAMPLER_MEAT_H +#define PG_SAMPLER_MEAT_H + + +//////////////////////////////////////////////////////////////////////////////// +// Gamma // +//////////////////////////////////////////////////////////////////////////////// + +template +inline T rpg_gamma(const T& shape_, const T& scale_, int trunc) { + PolyaGamma pg(trunc); + + // Assume + unsigned int n = shape_.size(); + + // Initialize a similar result object + T result(n); + + for(int i=0; i < n; ++i){ + +#ifdef USE_R + if (i % 1000 == 0) R_CheckUserInterrupt(); +#endif + + if (shape_[i] != 0.0) + result[i] = pg.draw_sum_of_gammas(shape_[i], scale_[i]); + else + result[i] = 0.0; + } + + return result; +} + +template +inline T rpg_gamma_trunc(const T& shape_, const T& scale_) { + // Set with default value to align function calls + return rpg_gamma(shape_, scale_, 1000); +} + +//////////////////////////////////////////////////////////////////////////////// +// Devroye // +//////////////////////////////////////////////////////////////////////////////// + +template +inline T rpg_devroye(const T& shape_, const T& scale_) { + PolyaGamma pg(1); + + // Assume + unsigned int n = shape_.size(); + + // Initialize a similar result object + T result(n); + + for(unsigned int i = 0; i < n; ++i){ + int shape = static_cast(shape_[i]); + if (shape != 0) + result[i] = pg.draw(shape, scale_[i]); + else + result[i] = 0.0; + } + + return result; +} + +//////////////////////////////////////////////////////////////////////////////// +// Alt // +//////////////////////////////////////////////////////////////////////////////// + + +// template +// inline T rpg_alt(const T& shape_, const T& scale_) { +// PolyaGammaApproxAlt pg; +// +// // Assume +// unsigned int n = shape_.size(); +// +// // Initialize a similar result object +// T result(n); +// +// for(unsigned int i = 0; i < n; ++i){ +// if (shape_[i]!=0) +// result[i] = pg.draw(shape_[i], scale_[i]); +// else +// result[i] = 0.0; +// } +// +// return result; +// } + +//////////////////////////////////////////////////////////////////////////////// +// Saddlepoint // +//////////////////////////////////////////////////////////////////////////////// + +template +inline T rpg_sp(const T& shape_, const T& scale_) { + PolyaGammaApproxSP pg_sp; + + // Assume + unsigned int n = shape_.size(); + + // Initialize a similar result object + T result(n); + + // Establish a double to handle return from pg_sp + double d = 0.0; + + for(unsigned int i = 0; i < n; ++i){ + if (shape_[i] != 0) { + // x[i] refers to an initialized vector with value 0. + pg_sp.draw(d, shape_[i], scale_[i]); + // Draw returns the number of iterations. + // Silent reference return of the double. + result[i] = d; + } else { + result[i] = 0.0; + } + } + + return result; +} + + +//////////////////////////////////////////////////////////////////////////////// +// Normal approximation // +//////////////////////////////////////////////////////////////////////////////// + +template +inline T rpg_normal(const T& shape_, const T& scale_) { + // Establish a Polya-Gamma class + PolyaGamma pg; + + // Assume + unsigned int n = shape_.size(); + + // Initialize a similar result object + T result(n); + + for(unsigned int i = 0; i < n; ++i) { + // Think shape is equivalent ot sample size + double m = pg.pg_m1(shape_[i], scale_[i]); + double v = pg.pg_m2(shape_[i], scale_[i]) - m * m; + result[i] = R::rnorm(m, sqrt(v)); + } + + return result; +} + +//////////////////////////////////////////////////////////////////////////////// +// Decider functions // +//////////////////////////////////////////////////////////////////////////////// + +inline double rpg_scalar_hybrid(double shape, double scale) { + + // Establish a Polya-Gamma class + PolyaGamma pg; + // Setup a Saddlepoint sampler + PolyaGammaApproxSP pg_sp; + + double result = 0.0; + + if (shape > 170) { + double m = pg.pg_m1(shape, scale); + double v = pg.pg_m2(shape, scale) - m*m; + result = R::rnorm(m, sqrt(v)); + } + else if (shape > 13) { + // x[i] refers to an initialized vector with value 0. + pg_sp.draw(result, shape, scale); + // Draw returns the number of iterations. + // Silent reference return of the double. + } + else if (shape == 1 || shape == 2) { + result = pg.draw(static_cast(shape), scale); + } + // Disabled. + // else if (shape > 1) { + // result = pg_alt.draw(shape, scale); + // } + else if (shape > 0) { + result = pg.draw_sum_of_gammas(shape, scale); + } + else { + result = 0.0; + } + + return result; +} + +inline arma::vec rpg_vector_hybrid( + unsigned int n, double shape, double scale) { + + // Establish a Polya-Gamma class + PolyaGamma pg; + // Setup a Saddlepoint sampler + PolyaGammaApproxSP pg_sp; + + arma::vec result(n); + + // Establish a double to handle return from pg_sp + double d = 0.0; + + for(unsigned int i = 0; i < n; ++i) { + if (shape > 170) { + double m = pg.pg_m1(shape, scale); + double v = pg.pg_m2(shape, scale) - m*m; + result[i] = R::rnorm(m, sqrt(v)); + } + else if (shape > 13) { + // x[i] refers to an initialized vector with value 0. + pg_sp.draw(d, shape, scale); + // Draw returns the number of iterations. + // Slient reference return of the double. + result[i] = d; + } + else if (shape == 1 || shape == 2) { + result[i] = pg.draw(static_cast(shape), scale); + } + // Disabled. + // else if (shape > 1) { + // result = pg_alt.draw(shape, scale); + // } + else if (shape > 0) { + result[i] = pg.draw_sum_of_gammas(shape, scale); + } + else { + result[i] = 0.0; + } + + } + + return result; +} + +// Assume shape_ and scale_ are equivalent +// Create a generic template update +template +inline T rpg_hybrid(const T& shape_, const T& scale_) { + + // Establish a Polya-Gamma class + PolyaGamma pg; + + // Setup a Saddlepoint sampler + PolyaGammaApproxSP pg_sp; + + // Assume + unsigned int n = shape_.size(); + + // Initialize a similar result object + T result(n); + + // Establish a double to handle return from pg_sp + double d = 0.0; + + for(unsigned int i = 0; i < n; ++i) { + double shape = shape_[i]; + double scale = scale_[i]; + + // Think shape is equivalent ot sample size + if (shape > 170) { + double m = pg.pg_m1(shape, scale); + double v = pg.pg_m2(shape, scale) - m * m; + result[i] = R::rnorm(m, sqrt(v)); + } + else if (shape > 13) { + // x[i] refers to an initialized vector with value 0. + pg_sp.draw(d, shape, scale); + // Draw returns the number of iterations. + // Silent reference return of the double. + result[i] = d; + } + else if (shape == 1 || shape == 2) { + result[i] = pg.draw(static_cast(shape), scale); + } + // Disabled. + // else if (shape > 1) { + // result = pg_alt.draw(shape, scale); + // } + else if (shape > 0) { + result[i] = pg.draw_sum_of_gammas(shape, scale); + } + else { + result[i] = 0.0; + } + + } + + return result; +} + +#endif diff --git a/inst/include/pg_truncated_gamma_bones.h b/inst/include/pg_truncated_gamma_bones.h new file mode 100644 index 0000000..60ff78c --- /dev/null +++ b/inst/include/pg_truncated_gamma_bones.h @@ -0,0 +1,33 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +// See Philippe, *Simulation of right and left truncated gamma +// distributions by mixtures* and Dagpunar *Sampling of variates from +// a truncated gamma distribution* + +#ifndef PG_TRUNCATED_GAMMA_BONES_H +#define PG_TRUNCATED_GAMMA_BONES_H + +inline double right_tgamma_reject(double shape, double rate); +inline double omega_k(int k, double a, double b); +inline double right_tgamma_beta(double shape, double rate); +inline double rtgamma_rate(double shape, double rate, double right_t); +inline double ltgamma(double shape, double rate, double trunc); + +#endif diff --git a/inst/include/pg_truncated_gamma_meat.h b/inst/include/pg_truncated_gamma_meat.h new file mode 100644 index 0000000..1343606 --- /dev/null +++ b/inst/include/pg_truncated_gamma_meat.h @@ -0,0 +1,126 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_TRUNCATED_GAMMA_MEAT_H +#define PG_TRUNCATED_GAMMA_MEAT_H + +// Truncatation at t = 1. +inline double right_tgamma_reject(double shape, double rate) +{ + double x = 2.0; + while (x > 1.0) + x = gamma_rate(shape, rate); + return x; +} + +double omega_k(int k, double a, double b) +{ + double log_coef = -b + (a+k-1) * log(b) - ln_gamma(a+k) - p_gamma_rate(1.0, a, b, true); + return exp(log_coef); +} + +// Truncation at t = 1. +double right_tgamma_beta(double shape, double rate) +{ + double a = shape; + double b = rate; + + double u = r_unif(); + + int k = 1; + double cdf = omega_k(1, a, b); + while (u > cdf) { + cdf += omega_k(++k, a, b); + if (k % 100000 == 0) { + #ifndef USE_R + fprintf(stderr, "right_tgamma_beta (itr k=%i): a=%g, b=%g, u=%g, cdf=%g\n", k, a, b, u, cdf); + #else + Rprintf("right_tgamma_beta (itr k=%i): a=%g, b=%g, u=%g, cdf=%g\n", k, a, b, u, cdf); + R_CheckUserInterrupt(); + #endif + } + } + + return R::beta(a, k); +} + +double rtgamma_rate(double shape, double rate, double right_t) +{ + // x \sim (a,b,t) + // ty = x + // y \sim (a, bt, 1); + double a = shape; + double b = rate * right_t; + + double p = p_gamma_rate(1.0, a, b, false); + double y = 0.0; + if (p > 0.95) + y = right_tgamma_reject(a, b); + else + y = right_tgamma_beta(a,b); + + double x = right_t * y; + return x; +} + +double ltgamma(double shape, double rate, double trunc) +{ + double a = shape; + double b = rate * trunc; + + if (trunc <=0) { + #ifndef USE_R + fprintf(stderr, "ltgamma: trunc = %g < 0\n", trunc); + #else + Rprintf("ltgamma: trunc = %g < 0\n", trunc); + #endif + return 0; + } + if (shape < 1) { + #ifndef USE_R + fprintf(stderr, "ltgamma: shape = %g < 1\n", shape); + #else + Rprintf("ltgamma: shape = %g < 1\n", shape); + #endif + return 0; + } + + if (shape ==1) return expon_rate(1.) / rate + trunc; + + double d1 = b-a; + double d3 = a-1; + double c0 = 0.5 * (d1 + sqrt(d1*d1 + 4 * b)) / b; + + double x = 0.0; + bool accept = false; + + while (!accept) { + x = b + expon_rate(1.) / c0; + double u = r_unif(); + + double l_rho = d3 * log(x) - x * (1-c0); + double l_M = d3 * log(d3 / (1-c0)) - d3; + + accept = log(u) <= (l_rho - l_M); + } + + return trunc * (x/b); +} + +#endif diff --git a/inst/include/pg_truncated_norm_bones.h b/inst/include/pg_truncated_norm_bones.h new file mode 100644 index 0000000..d22f2f2 --- /dev/null +++ b/inst/include/pg_truncated_norm_bones.h @@ -0,0 +1,50 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +// Method from Christian Robert, *Simulation of truncated normal variables* + +// Code modified from + +#ifndef PG_TRUNCATED_NORM_BONES_H +#define PG_TRUNCATED_NORM_BONES_H + +// Throw runtime exception or print message and return. +#ifndef TREOR +#ifndef NTHROW +#define TREOR(MESS, VAL) throw std::runtime_error(MESS); +#else +#ifndef USE_R +#define TREOR(MESS, VAL) {fprintf(stderr, MESS); return VAL;} +#else +#define TREOR(MESS, VAL) {Rprintf(MESS); return VAL;} +#endif +#endif +#endif + +inline double flat(double left, double right); +inline double texpon_rate(double left, double rate); +inline double texpon_rate(double left, double right, double rate); +inline double alphastar(double left); +inline double lowerbound(double left); +inline double tnorm(double left); +inline double tnorm(double left, double mu, double sd); +inline double tnorm(double left, double right, double mu, double sd); +inline double rtinvchi2(double scale, double trunc); + +#endif diff --git a/inst/include/pg_truncated_norm_meat.h b/inst/include/pg_truncated_norm_meat.h new file mode 100644 index 0000000..eab08d5 --- /dev/null +++ b/inst/include/pg_truncated_norm_meat.h @@ -0,0 +1,275 @@ +// This header file contains work done by +// (C) Nicholas Polson, James Scott, Jesse Windle, 2012-2019 +// within the BayesLogit R package. This work is licensed under GPL v3 or greater. + +// Modifications for header-only inclusion was done by +// (C) James Balamuta 2019 - 2020 + +// pg is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation, either version 3 of the License, or (at your option) any later +// version. + +// pg is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License along with +// pg. If not, see . + +#ifndef PG_TRUNCATED_NORM_MEAT_H +#define PG_TRUNCATED_NORM_MEAT_H + +inline double flat(double a, double b) { + return a + r_unif() * (b-a); +} + +inline double texpon_rate(double left, double rate){ + if (rate < 0) TREOR("texpon_rate: rate < 0, return 0\n", 0.0); + // return left - log(r_unif()) / rate; + return expon_rate(rate) + left; +} + +inline double texpon_rate(double left, double right, double rate) +{ + if (left == right) return left; + if (left > right) TREOR("texpon_rate: left > right, return 0.\n", 0.0); + if (rate < 0) TREOR("texpon_rate: rate < 0, return 0\n", 0.0); + + double b = 1 - exp(rate * (left - right)); + double y = 1 - b * r_unif(); + return left - log(y) / rate; +} + +inline double alphastar(double left) +{ + return 0.5 * (left + sqrt(left*left + 4)); +} + +inline double lowerbound(double left) +{ + double astar = alphastar(left); + double lbound = left + exp(0.5 * + 0.5 * left * (left - astar)) / astar; + return lbound; +} + +inline double tnorm(double left) +{ + double rho, ppsl; + int count = 1; + + if (left < 0) { // Accept/Reject Normal + while (true) { + ppsl = r_norm(0.0, 1.0); + if (ppsl > left) return ppsl; + check_R_interupt(count++); + #ifndef NDEBUG + if (count > RCHECK * 1000) { + #ifndef USE_R + fprintf(stderr, "left < 0; count: %i\n", count); + #else + Rprintf("left < 0; count: %i\n", count); + #endif + } + //p + #endif + } + } + else { // Accept/Reject Exponential + // return tnorm_tail(left); // Use Devroye. + double astar = alphastar(left); + while (true) { + ppsl = texpon_rate(left, astar); + rho = exp( -0.5 * (ppsl - astar) * (ppsl - astar) ); + if (r_unif() < rho) return ppsl; + check_R_interupt(count++); + #ifndef NDEBUG + if (count > RCHECK * 1000) { + #ifndef USE_R + fprintf(stderr, "left > 0; count: %i\n", count); + #else + Rprintf("left > 0; count: %i\n", count); + #endif + } + #endif + } + } + +} + +inline double tnorm(double left, double right) +{ + // The most difficult part of this algorithm is figuring out all the + // various cases. An outline is summarized in the Appendix. + + // Check input + #ifdef USE_R + if (ISNAN(right) || ISNAN(left)) + #else + if (std::isnan(right) || std::isnan(left)) + #endif + { + #ifndef USE_R + fprintf(stderr, "Warning: nan sent to tnorm: left=%g, right=%g\n", left, right); + #else + Rprintf("Warning: nan sent to tnorm: left=%g, right=%g\n", left, right); + #endif + TREOR("tnorm: parameter problem.\n", 0.5 * (left + right)); + // throw std::runtime_error("tnorm: parameter problem.\n"); + } + + if (right < left) { + #ifndef USE_R + fprintf(stderr, "Warning: left: %g, right:%g.\n", left, right); + #else + Rprintf("Warning: left: %g, right:%g.\n", left, right); + #endif + TREOR("tnorm: parameter problem.\n", 0.5 * (left + right)); + } + + double rho, ppsl; + int count = 1; + + if (left >= 0) { + double lbound = lowerbound(left); + if (right > lbound) { // Truncated Exponential. + double astar = alphastar(left); + while (true) { + ppsl = texpon_rate(left, right, astar); + rho = exp(-0.5*(ppsl - astar)*(ppsl-astar)); + if (r_unif() < rho) return ppsl; + if (count > RCHECK * 10) { + #ifndef USE_R + fprintf(stderr, "left >= 0, right > lbound; count: %i\n", count); + #else + Rprintf("left >= 0, right > lbound; count: %i\n", count); + #endif + } + // if (ppsl < right) return ppsl; + } + } + else { + while (true) { + ppsl = flat(left, right); + rho = exp(0.5 * (left*left - ppsl*ppsl)); + if (r_unif() < rho) return ppsl; + check_R_interupt(count++); + #ifndef NDEBUG + if (count > RCHECK * 10) { + #ifndef USE_R + fprintf(stderr, "left >= 0, right <= lbound; count: %i\n", count); + #else + Rprintf("left >= 0, right <= lbound; count: %i\n", count); + #endif + } + #endif + } + } + } + else if (right >= 0) { + if ( (right - left) < SQRT2PI ){ + while (true) { + ppsl = flat(left, right); + rho = exp(-0.5 * ppsl * ppsl); + if (r_unif() < rho) return ppsl; + check_R_interupt(count++); + #ifndef NDEBUG + if (count > RCHECK * 10) { + #ifndef USE_R + fprintf(stderr, "First, left < 0, right >= 0, count: %i\n", count); + #else + Rprintf("First, left < 0, right >= 0, count: %i\n", count); + #endif + } + #endif + } + } + else{ + while (true) { + ppsl = r_norm(0., 1.); + if (left < ppsl && ppsl < right) return ppsl; + check_R_interupt(count++); + #ifndef NDEBUG + if (count > RCHECK * 10) { + #ifndef USE_R + fprintf(stderr, "Second, left < 0, right > 0, count: %i\n", count); + #else + Rprintf("Second, left < 0, right > 0, count: %i\n", count); + #endif + } + #endif + } + } + } + else { + return -1. * tnorm(-1.0 * right, -1.0 * left); + } +} + +inline double tnorm(double left, double mu, double sd) +{ + double newleft = (left - mu) / sd; + return mu + tnorm(newleft) * sd; +} + +inline double tnorm(double left, double right, double mu, double sd) +{ + if (left==right) return left; + + double newleft = (left - mu) / sd; + double newright = (right - mu) / sd; + + // I want to check this here as well so we can see what the input was. + // It may be more elegant to try and catch tdraw. + if (newright < newleft) { + #ifndef USE_R + fprintf(stderr, "left, right, mu, sd: %g, %g, %g, %g \n", left, right, mu, sd); + fprintf(stderr, "nleft, nright: %g, %g\n", newleft, newright); + #else + Rprintf("left, right, mu, sd: %g, %g, %g, %g \n", left, right, mu, sd); + Rprintf("nleft, nright: %g, %g\n", newleft, newright); + #endif + TREOR("tnorm: parameter problem.\n", 0.5 * (left + right)); + } + + double tdraw = tnorm(newleft, newright); + double draw = mu + tdraw * sd; + + // It may be the case that there is some numerical error and that the draw + // ends up out of bounds. + if (draw < left || draw > right){ + #ifndef USE_R + fprintf(stderr, "Error in tnorm: draw not in bounds.\n"); + fprintf(stderr, "left, right, mu, sd: %g, %g, %g, %g\n", left, right, mu, sd); + fprintf(stderr, "nleft, nright, tdraw, draw: %g, %g, %g, %g\n", newleft, newright, tdraw, draw); + #else + Rprintf("Error in tnorm: draw not in bounds.\n"); + Rprintf("left, right, mu, sd: %g, %g, %g, %g\n", left, right, mu, sd); + Rprintf("nleft, nright, tdraw, draw: %g, %g, %g, %g\n", newleft, newright, tdraw, draw); + #endif + TREOR("Aborting and returning average of left and right.\n", 0.5 * (left + right)); + } + + return draw; +} + +inline double rtinvchi2(double scale, double trunc) +{ + double R = trunc / scale; + // double X = 0.0; + // // I need to consider using a different truncated normal sampler. + // double E1 = r.expon_rate(1.0); double E2 = r.expon_rate(1.0); + // while ( (E1*E1) > (2 * E2 / R)) { + // // printf("E %g %g %g %g\n", E1, E2, E1*E1, 2*E2/R); + // E1 = r.expon_rate(1.0); E2 = r.expon_rate(1.0); + // } + // // printf("E %g %g \n", E1, E2); + // X = 1 + E1 * R; + // X = R / (X * X); + // X = scale * X; + double E = tnorm(1/sqrt(R)); + double X = scale / (E*E); + return X; +} + +#endif diff --git a/man/pg-package.Rd b/man/pg-package.Rd new file mode 100644 index 0000000..92adb23 --- /dev/null +++ b/man/pg-package.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pg-package.R +\docType{package} +\name{pg-package} +\alias{pg} +\alias{pg-package} +\title{pg: Polya-Gamma Distribution Sampler} +\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'. +} +\author{ +\strong{Maintainer}: James Balamuta \email{balamut2@illinois.edu} [copyright holder] + +} +\keyword{internal} diff --git a/man/rpg-sampler.Rd b/man/rpg-sampler.Rd new file mode 100644 index 0000000..754a111 --- /dev/null +++ b/man/rpg-sampler.Rd @@ -0,0 +1,92 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{rpg_scalar} +\alias{rpg_scalar} +\alias{rpg_vector} +\alias{rpg_hybrid} +\alias{rpg_gamma} +\alias{rpg_devroye} +\alias{rpg_sp} +\alias{rpg_normal} +\title{Sample from the Polya-Gamma distribution PG(h, z)} +\usage{ +rpg_scalar(h, z) + +rpg_vector(n, h, z) + +rpg_hybrid(h, z) + +rpg_gamma(h, z, trunc = 1000L) + +rpg_devroye(h, z) + +rpg_sp(h, z) + +rpg_normal(h, z) +} +\arguments{ +\item{h}{\code{integer} values corresponding to the "shape" parameter.} + +\item{z}{\code{numeric} values corresponding to the "scale" parameter.} + +\item{n}{The number of samples to taken from a PG(h, z). Used only by +the vector sampler.} + +\item{trunc}{Truncation cut-off. Only used by the gamma sampler.} +} +\value{ +A single \code{numeric} value. +} +\description{ +Chooses the most efficient implemented method to sample from a Polya-Gamma +distribution. Details on algorithm selection presented below. +} +\details{ +The following sampling cases are enabled: +\itemize{ +\item \code{h > 170}: Normal approximation method +\item \code{h > 13}: Saddlepoint approximation method +\item \code{h = 1} or \code{h = 2}: Devroye method +\item \code{h > 0}: Sum of Gammas method. +\item \code{h < 0}: Result is automatically set to zero. +} +} +\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) +} diff --git a/man/theoretical-pg.Rd b/man/theoretical-pg.Rd new file mode 100644 index 0000000..f784255 --- /dev/null +++ b/man/theoretical-pg.Rd @@ -0,0 +1,48 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pg_theoretical.R +\name{pg_mean} +\alias{pg_mean} +\alias{pg_var} +\title{Theoretical Polya-Gamma Distribution's Mean and Variance} +\usage{ +pg_mean(h, z) + +pg_var(h, z) +} +\arguments{ +\item{h}{A single \code{integer} value corresponding to the "shape" parameter.} + +\item{z}{A single \code{numeric} value corresponding to the "scale" parameter.} +} +\value{ +Either the theoretical mean or theoretical variance for a Polya-Gamma +distribution. +} +\description{ +Compute the theoretical mean and variance for a Polya-Gamma variable. +} +\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)) +} diff --git a/pg.Rproj b/pg.Rproj new file mode 100644 index 0000000..88ff2b5 --- /dev/null +++ b/pg.Rproj @@ -0,0 +1,21 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: knitr +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source +PackageRoxygenize: rd,collate,namespace diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..22034c4 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,3 @@ +*.o +*.so +*.dll diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..678ba52 --- /dev/null +++ b/src/Makevars @@ -0,0 +1 @@ +PKG_CXXFLAGS=-I../inst/include/ diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 0000000..678ba52 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1 @@ +PKG_CXXFLAGS=-I../inst/include/ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..ae4b5bb --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,129 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include "../inst/include/pg.h" +#include +#include + +using namespace Rcpp; + +#ifdef RCPP_USE_GLOBAL_ROSTREAM +Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +#endif + +// rpg_scalar +double rpg_scalar(const double h, const double z); +RcppExport SEXP _pg_rpg_scalar(SEXP hSEXP, SEXP zSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const double >::type h(hSEXP); + Rcpp::traits::input_parameter< const double >::type z(zSEXP); + rcpp_result_gen = Rcpp::wrap(rpg_scalar(h, z)); + return rcpp_result_gen; +END_RCPP +} +// rpg_vector +arma::vec rpg_vector(unsigned int n, const double h, const double z); +RcppExport SEXP _pg_rpg_vector(SEXP nSEXP, SEXP hSEXP, SEXP zSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP); + Rcpp::traits::input_parameter< const double >::type h(hSEXP); + Rcpp::traits::input_parameter< const double >::type z(zSEXP); + rcpp_result_gen = Rcpp::wrap(rpg_vector(n, h, z)); + return rcpp_result_gen; +END_RCPP +} +// rpg_hybrid +arma::vec rpg_hybrid(const arma::vec& h, const arma::vec& z); +RcppExport SEXP _pg_rpg_hybrid(SEXP hSEXP, SEXP zSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type h(hSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + rcpp_result_gen = Rcpp::wrap(rpg_hybrid(h, z)); + return rcpp_result_gen; +END_RCPP +} +// rpg_scalar_loop +arma::vec rpg_scalar_loop(const arma::vec& h, const arma::vec& z); +RcppExport SEXP _pg_rpg_scalar_loop(SEXP hSEXP, SEXP zSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type h(hSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + rcpp_result_gen = Rcpp::wrap(rpg_scalar_loop(h, z)); + return rcpp_result_gen; +END_RCPP +} +// rpg_gamma +arma::vec rpg_gamma(const arma::vec& h, const arma::vec& z, int trunc); +RcppExport SEXP _pg_rpg_gamma(SEXP hSEXP, SEXP zSEXP, SEXP truncSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type h(hSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + Rcpp::traits::input_parameter< int >::type trunc(truncSEXP); + rcpp_result_gen = Rcpp::wrap(rpg_gamma(h, z, trunc)); + return rcpp_result_gen; +END_RCPP +} +// rpg_devroye +arma::vec rpg_devroye(const arma::vec& h, const arma::vec& z); +RcppExport SEXP _pg_rpg_devroye(SEXP hSEXP, SEXP zSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type h(hSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + rcpp_result_gen = Rcpp::wrap(rpg_devroye(h, z)); + return rcpp_result_gen; +END_RCPP +} +// rpg_sp +arma::vec rpg_sp(const arma::vec& h, const arma::vec& z); +RcppExport SEXP _pg_rpg_sp(SEXP hSEXP, SEXP zSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type h(hSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + rcpp_result_gen = Rcpp::wrap(rpg_sp(h, z)); + return rcpp_result_gen; +END_RCPP +} +// rpg_normal +arma::vec rpg_normal(const arma::vec& h, const arma::vec& z); +RcppExport SEXP _pg_rpg_normal(SEXP hSEXP, SEXP zSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type h(hSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + rcpp_result_gen = Rcpp::wrap(rpg_normal(h, z)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_pg_rpg_scalar", (DL_FUNC) &_pg_rpg_scalar, 2}, + {"_pg_rpg_vector", (DL_FUNC) &_pg_rpg_vector, 3}, + {"_pg_rpg_hybrid", (DL_FUNC) &_pg_rpg_hybrid, 2}, + {"_pg_rpg_scalar_loop", (DL_FUNC) &_pg_rpg_scalar_loop, 2}, + {"_pg_rpg_gamma", (DL_FUNC) &_pg_rpg_gamma, 3}, + {"_pg_rpg_devroye", (DL_FUNC) &_pg_rpg_devroye, 2}, + {"_pg_rpg_sp", (DL_FUNC) &_pg_rpg_sp, 2}, + {"_pg_rpg_normal", (DL_FUNC) &_pg_rpg_normal, 2}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_pg(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/rpg.cpp b/src/rpg.cpp new file mode 100644 index 0000000..b5330b0 --- /dev/null +++ b/src/rpg.cpp @@ -0,0 +1,125 @@ +#include "pg.h" + +//' 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) +// [[Rcpp::export]] +double rpg_scalar(const double h, const double z) { + return pg::rpg_scalar_hybrid(h, z); +} + +//' @export +//' @rdname rpg-sampler +// [[Rcpp::export]] +arma::vec rpg_vector(unsigned int n, const double h, const double z) { + return pg::rpg_vector_hybrid(n, h, z); +} + +//' @export +//' @rdname rpg-sampler +// [[Rcpp::export]] +arma::vec rpg_hybrid(const arma::vec& h, const arma::vec& z) { + return pg::rpg_hybrid(h, z); +} + + +// Local function only available inside of the package namespace. +// [[Rcpp::export]] +arma::vec rpg_scalar_loop(const arma::vec& h, const arma::vec& z) { + unsigned int n = h.size(); + arma::vec result(n); + for(unsigned int i = 0; i < h.size(); ++i) { + result[i] = pg::rpg_scalar_hybrid(h[i], z[i]); + } + + return result; +} + +//' @export +//' @rdname rpg-sampler +// [[Rcpp::export]] +arma::vec rpg_gamma(const arma::vec& h, const arma::vec& z, int trunc = 1000) { + return pg::rpg_gamma(h, z, trunc); +} + +//' @export +//' @rdname rpg-sampler +// [[Rcpp::export]] +arma::vec rpg_devroye(const arma::vec& h, const arma::vec& z) { + return pg::rpg_devroye(h, z); +} + +//' @export +//' @rdname rpg-sampler +// [[Rcpp::export]] +arma::vec rpg_sp(const arma::vec& h, const arma::vec& z) { + return pg::rpg_sp(h, z); +} + +//' @export +//' @rdname rpg-sampler +// [[Rcpp::export]] +arma::vec rpg_normal(const arma::vec& h, const arma::vec& z) { + return pg::rpg_normal(h, z); +} diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..2358867 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(pg) + +test_check("pg") diff --git a/tests/testthat/test-saddlepoint.R b/tests/testthat/test-saddlepoint.R new file mode 100644 index 0000000..a3b5fc9 --- /dev/null +++ b/tests/testthat/test-saddlepoint.R @@ -0,0 +1,49 @@ +test_that("Verify repeated sampling under same parameters matches", { + + h = rep_len(c(1, 2, 3), 100) + z = rep_len(c(4, 5, 6), 100) + + set.seed(5123) + x1 = pg::rpg_sp(h, z) + set.seed(5123) + x2 = pg::rpg_sp(h, z) + + expect_equal(x1, x2, info = "Samples should be equivalent") +}) + +test_that("Ensure setting a seed results in different values being generated", { + h = rep_len(c(1, 2, 3), 100) + z = rep_len(c(4, 5, 6), 100) + + set.seed(9143) + x1 = pg::rpg_sp(h, z) + set.seed(5123) + x2 = pg::rpg_sp(h, z) + + expect_false(isTRUE(all.equal(x1, x2)), info = "Samples should differ under seeds") +}) + +test_that("Ensure estimation is reasonably matching with theoretical mean", { + h = 100 + z = 4 + + h_vec = rep_len(h, 1e6) + z_vec = rep_len(z, 1e6) + + set.seed(5503) + x = pg::rpg_sp(h_vec, z_vec) + + mean_empirical_recovery = mean(x) + mean_theoretical_recovery = pg::pg_mean(h, z) + + var_empirical_recovery = c(var(x)) + var_theoretical_recovery = pg::pg_var(h, z) + + expect_equal(mean_empirical_recovery, mean_theoretical_recovery, + tolerance = 0.01, + info = "Sample mean matches") + + expect_equal(var_empirical_recovery, var_theoretical_recovery, + tolerance = 0.01, + info = "Sample variance matches") +}) diff --git a/tests/testthat/test-sample-control.R b/tests/testthat/test-sample-control.R new file mode 100644 index 0000000..afc9779 --- /dev/null +++ b/tests/testthat/test-sample-control.R @@ -0,0 +1,46 @@ +test_that("Ensure scalar draws match under given seed", { + set.seed(1512) + x1 = pg::rpg_scalar(1, 0.5) + set.seed(1512) + x2 = pg::rpg_scalar(1, 0.5) + + expect_equal(x1, x2, info = "Both samples should be equivalent") +}) + +test_that("Check scalar draws match reused C++ class draws", { + set.seed(1512) + + h = c(1, 25, 35, 50, 100, 171) + z = c(0.5, 4, 5, 6, 7, 8) + + x1 = rep(NA, length(h)) + for(i in seq_along(x1)) { + x1[i] = pg::rpg_scalar(h[i], z[i]) + } + # Cast to matrix + x1 = matrix(x1) + + set.seed(1512) + x2 = pg::rpg_hybrid(h, z) + + expect_equal(x1, x2, info = "Both samples should be equivalent") +}) + + +test_that("Verify repeated sampling under same parameters matches", { + set.seed(5123) + x1 = pg::rpg_vector(10, 1, 0.5) + set.seed(5123) + x2 = pg::rpg_vector(10, 1, 0.5) + + expect_equal(x1, x2, info = "Both vectorized samples should be equivalent") +}) + +test_that("Ensure setting a seed results in different values being generated", { + set.seed(9321) + x1 = pg::rpg_vector(10, 1, 0.5) + set.seed(5123) + x2 = pg::rpg_vector(10, 1, 0.5) + + expect_false(isTRUE(all.equal(x1, x2)), info = "Samples should differ under seeds") +})