The pgdraw package

Description

This package contains a function to generates random samples from the Polya-Gamma distribution using an implementation of the algorithm described in J. Windle's PhD thesis. A frequent application of this distribution is in Bayesian analysis of logistic regression models.

Details

The underlying implementation is in C.

For usage, see the examples in pgdraw and pgdraw.moments.

Note

To cite this package please reference:

Makalic, E. & Schmidt, D. F. High-Dimensional Bayesian Regularised Regression with the BayesReg Package arXiv:1611.06649 [stat.CO], 2016 https://arxiv.org/pdf/1611.06649.pdf

A MATLAB-compatible implementation of the sampler in this package can be obtained from:

http://dschmidt.org/?page_id=189

Author(s)

Daniel Schmidt daniel.schmidt@monash.edu

Faculty of Information Technology, Monash University, Australia

Enes Makalic emakalic@unimelb.edu.au

Centre for Epidemiology and Biostatistics, The University of Melbourne, Australia

References

Jesse Bennett Windle Forecasting High-Dimensional, Time-Varying Variance-Covariance Matrices with High-Frequency Data and Sampling Polya-Gamma Random Variates for Posterior Distributions Derived from Logistic Likelihoods PhD Thesis, 2013

Bayesian Inference for Logistic Models Using Polya-Gamma Latent Variables Nicholas G. Polson, James G. Scott and Jesse Windle Journal of the American Statistical Association Vol. 108, No. 504, pp. 1339–1349, 2013

Chung, Y.: Simulation of truncated gamma variables Korean Journal of Computational & Applied Mathematics, 1998, 5, 601-610

See Also

pgdraw, pgdraw.moments

Generate random samples from the Polya-Gamma distribution, PG(b,c)

Description

Generate random samples from the Polya-Gamma distribution

Usage

pgdraw(b, c)

Arguments

Value

A vector of samples from the Polya-Gamma distribution, one for each entry of c

Details

This code generates random variates from the Polya-Gamma distribution with desired 'b' and 'c' parameters. The underlying code is written in C and is an implementation of the algorithm described in J. Windle's PhD thesis.

The main application of the Polya-Gamma distribution is in Bayesian analysis as it allows for a data augmentation (via a scale mixture of normals) approach for representation of the logistic regression likelihood (see Example 2 below).

Note

To cite this package please reference:

Makalic, E. & Schmidt, D. F. High-Dimensional Bayesian Regularised Regression with the BayesReg Package arXiv:1611.06649 [stat.CO], 2016 https://arxiv.org/pdf/1611.06649.pdf

A MATLAB-compatible implementation of the sampler in this package can be obtained from:

http://dschmidt.org/?page_id=189

References

Jesse Bennett Windle Forecasting High-Dimensional, Time-Varying Variance-Covariance Matrices with High-Frequency Data and Sampling Polya-Gamma Random Variates for Posterior Distributions Derived from Logistic Likelihoods, PhD Thesis, 2013

Bayesian Inference for Logistic Models Using Polya-Gamma Latent Variables Nicholas G. Polson, James G. Scott and Jesse Windle, Journal of the American Statistical Association Vol. 108, No. 504, pp. 1339–1349, 2013

Chung, Y.: Simulation of truncated gamma variables, Korean Journal of Computational & Applied Mathematics, 1998, 5, 601-610

See Also

pgdraw.moments

Examples

# -----------------------------------------------------------------
# Example 1: Simulated vs exact moments
u = matrix(1,1e6,1)
x = pgdraw(1,0.5*u)
mean(x)
var(x)
pgdraw.moments(1,0.5)

x = pgdraw(2,2*u)
mean(x)
var(x)
pgdraw.moments(2,2)


# -----------------------------------------------------------------
# Example 2: Simple logistic regression
#   Sample from the following Bayesian hierarchy:
#    y_i ~ Be(1/(1+exp(-b)))
#    b   ~ uniform on R (improper)
#
#   which is equivalent to
#    y_i - 1/2 ~ N(b, 1/omega2_i)
#    omega2_i  ~ PG(1,0)
#    b         ~ uniform on R
#
sample_simple_logreg <- function(y, nsamples)
{
  n = length(y)
  omega2 = matrix(1,n,1)   # Polya-Gamma latent variables
  beta   = matrix(0,nsamples,1)
  
  for (i in 1:nsamples)
  {
    # Sample 'beta'
    s = sum(omega2)
    m = sum(y-1/2)/s
    beta[i] = rnorm(1, m, sqrt(1/s))
    # Sample P-G L.Vs
    omega2 = pgdraw(1, matrix(1,n,1)*beta[i])
  }
  return(beta)
}

# 3 heads, 7 tails; ML estimate of p = 3/10 = 0.3
y = c(1,1,1,0,0,0,0,0,0,0)  

# Sample
b = sample_simple_logreg(y, 1e4)
hist(x=b)

# one way of estimating of 'p' from posterior samples
1/(1+exp(-mean(b)))         

Compute exact first and second moments for the Polya-Gamma distribution, PG(b, c)

Description

Compute exact first and second moments for the Polya-Gamma distribution

Usage

pgdraw.moments(b, c)

Arguments

Value

A list containing the mean and variance.

Details

This code computes the exact mean and variance of the Polya-Gamma distribution for the specified parameters.

References

Jesse Bennett Windle Forecasting High-Dimensional, Time-Varying Variance-Covariance Matrices with High-Frequency Data and Sampling Polya-Gamma Random Variates for Posterior Distributions Derived from Logistic Likelihoods PhD Thesis, 2013

Bayesian Inference for Logistic Models Using Polya-Gamma Latent Variables Nicholas G. Polson, James G. Scott and Jesse Windle Journal of the American Statistical Association Vol. 108, No. 504, pp. 1339–1349, 2013

See Also

pgdraw

Examples

# -----------------------------------------------------------------
# Example: Simulated vs exact moments

u = matrix(1,1e6,1)
x = pgdraw(1,0.5*u)
mean(x)
var(x)
pgdraw.moments(1,0.5)

x = pgdraw(2,2*u)
mean(x)
var(x)
pgdraw.moments(2,2)