-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.Rmd
167 lines (126 loc) · 5.14 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%"
)
```
# pg
<!-- badges: start -->
[](https://github.com/tmsalab/pg/actions/workflows/R-CMD-check.yaml)
[-brightgreen.svg?style=flat)](http://www.gnu.org/licenses/gpl-3.0.html)
<!-- badges: end -->
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 <pg.h>
// [[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://doi.org/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)
<https://arxiv.org/abs/1405.0506>. 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)