Translating Probability Density Functions: From R to BUGS and Back Again

The ability to implement statistical models in the BUGS language facilitates Bayesian inference by automating MCMC algorithms. Software packages that interpret the BUGS language include OpenBUGS, WinBUGS, and JAGS. R packages that link BUGS software to the R environment, including rjags and R2WinBUGS, are widely used in Bayesian analysis. Indeed, many packages in the Bayesian task view on CRAN (http://cran.r-project.org/web/views/Bayesian.html) depend on this integration. However, the R and BUGS languages use different representations of common probability density functions, creating a potential for errors to occur in the implementation or interpretation of analyses that use both languages. Here we review different parameterizations used by the R and BUGS languages, describe how to translate between the languages, and provide an R function, r2bugs.distributions, that transforms parameterizations from R to BUGS and back again.

Table 1: Summary of different parameterizations of common distributions used by R and BUGS. Note: For ease of reference, parameterizations follow the JAGS and R documentation; as a result, the table includes equivalent equations that appear different, either because JAGS and R use different names for the same parameter or because the equation has been rearranged. For example, the shape parameter of the Gamma distribution is \(r\) in the BUGS documentation and \(a\) in the R documentation. For the Binomial, Negative Binomial, and Gamma distributions, BUGS and R expect parameters in different order (the order of parameters matters since arguments are assigned based on position in BUGS and may be in R as well). R allows alternate parameterizations for the Negative Binomial and Gamma distributions, but these are not shown here. The variable \(x\) is implicit in all of the BUGS “Use” expressions. The Beta, Poisson, Exponential, and Uniform distributions have identical parameterizations in R and BUGS.

Distribution	Lang.	Parameterization	Use	Notes
Normal
	R	\(\frac{1}{\sqrt{2 \pi}\sigma}\exp\left(-\frac{\left(x - \mu\right)^2}{2 \sigma^2}\right)\)	dnorm(\(x\), mean =\(\mu\), sd =\(\sigma\))
	BUGS	\(\sqrt{\frac{\tau}{2\pi}}\exp\left(-\left(x-\mu\right)^2\tau\right)\)	dnorm(mean =\(\mu\), precision =\(\tau\))	\(\tau=\left(\frac{1}{\sigma}\right)^2\)
log-Normal
	R	\(\frac{1}{\sqrt{2 \pi} \sigma x} \exp\left(-\frac{\left(\textrm{log}\left(x\right) - \mu\right)^2}{\left(2 \sigma^2\right)}\right)\)	dlnorm(\(x\), mean =\(\mu\), sd =\(\sigma\))
	BUGS	\(\frac{\sqrt{\tau}}{x}\exp\left(\frac{-\tau\left(\textrm{log}\left(x\right)-\mu\right)^2}{2}\right)\)	dlnorm(mean =\(\mu\), precision =\(\tau\))	\(\tau=\left(\frac{1}{\sigma}\right)^2\)
Binomial				reverse parameter order
	R	\({n \choose x} p^{x}\left(1-p\right)^{n-x}\)	dbinom(\(x\), size =\(n\), prob =\(p\))
	BUGS	same	dbin(prob =\(p\), size =\(n\))
Negative Binomial				reverse parameter order
	R	\(\frac{\Gamma\left(x+n\right)}{\Gamma\left(n\right) x!} p^n \left(1-p\right)^x\)	dnbinom(\(x\), size =\(n\), prob =\(p\))	size (n) is continuous
	BUGS	\({x+r-1 \choose x}p^r\left(1-p\right)^x\)	dnegbin(prob =\(p\), size =\(r\))	size (r) is discrete
Weibull
	R	\(\frac{a}{b} (\frac{x}{b})^{a-1} \exp\left(- \left(\frac{x}{b}\right)^a\right)\)	dweibull(\(x\), shape =\(a\), scale =\(b\))
	BUGS	\(\nu\lambda x^{\nu - 1}\exp\left(-\lambda x^{\nu}\right)\)	dweib(shape =\(\nu\), lambda =\(\lambda\))	\(\lambda=\left(\frac{1}{b}\right)^a\)
Gamma				reverse parameter order
	R	\({\frac{r^a}{\Gamma(a)}} x^{a-1} \exp(-xr)\)	dgamma(\(x\), shape =\(a\), rate =\(r\))
	BUGS	\({\frac{\lambda^r x^{r-1}\exp(-\lambda x)}{\Gamma(r)}}\)	dgamma(shape =\(r\), lambda =\(\lambda\))

1 Probability density functions in R and BUGS

R and BUGS implement many of the same probability distribution functions, but they often parameterize the same distribution differently (Table 1). Although these probability distribution functions are clearly described in the documentation of their respective languages, we were unable to find a summary of these differences in one place. The motivation for this article is to document and clarify these differences. Our sources are the JAGS documentation (Plummer 2010) and the documentation of individual R functions.

2 A bilingual translation function

To support the automation of model specification in JAGS with priors computed and stored in R (LeBauer et al. 2013), we developed a function to translate parameterizations of common probability distributions from R to BUGS (and back again, by specifying direction = ’bugs2r’). Parameter transformations, parameter order, and differences in function names are documented in Table 1 and implemented in the R function r2bugs.distributions.

r2bugs.distributions <- function(priors, direction = 'r2bugs') {
  priors$distn  <- as.character(priors$distn)
  priors$parama <- as.numeric(priors$parama)
  priors$paramb <- as.numeric(priors$paramb)
  ## index dataframe according to distribution
  norm   <- priors$distn %in% c('norm', 'lnorm')    # these have same transform
  weib   <- grepl("weib", priors$distn)             # matches r and bugs version
  gamma  <- priors$distn == 'gamma'
  chsq   <- grepl("chisq", priors$distn)            # matches r and bugs version
  bin    <- priors$distn %in% c('binom', 'bin')     # matches r and bugs version
  nbin   <- priors$distn %in% c('nbinom', 'negbin') # matches r and bugs version
  
  ## Normal, log-Normal: Convert sd to precision
  exponent <- ifelse(direction == "r2bugs", -2, -0.5) 
  priors$paramb[norm] <-  priors$paramb[norm] ^ exponent
  
  ## Weibull
  if(direction == 'r2bugs'){
    ## Convert R parameter b to BUGS parameter lambda by l = (1/b)^a
    priors$paramb[weib] <- (1 / priors$paramb[weib]) ^ priors$parama[weib]
  } else if (direction == 'bugs2r') {
    ## Convert BUGS parameter lambda to BUGS parameter b by b = l^(-1/a)
    priors$paramb[weib] <-  priors$paramb[weib] ^ (- 1 / priors$parama[weib] ) 
  }
  
  ## Reverse parameter order for binomial and negative binomial
  priors[bin | nbin, c('parama', 'paramb')] <-
    priors[bin | nbin, c('paramb', 'parama')]
  
  ## Translate distribution names
  if(direction == "r2bugs"){
    priors$distn[weib] <- "weib"
    priors$distn[chsq] <- "chisqr"
    priors$distn[bin]  <- "bin"
    priors$distn[nbin] <- "negbin"
  } else if(direction == "bugs2r"){
    priors$distn[weib] <- "weibull"
    priors$distn[chsq] <- "chisq"
    priors$distn[bin]  <- "binom"
    priors$distn[nbin] <- "nbinom"
  }
  return(priors)
}

3 A simple example

As an example, we take the R-parameterized prior distribution \(X \sim \mathcal{N}(\mu=10,\sigma=2)\) and convert it to BUGS parameterization \(X \sim \mathcal{N}(\mu=10,\tau=1/4)\). We specify a model in JAGS that allows us to sample directly from a prior distribution. The function works for each of the distributions in Table 1. This particular example is the JAGS implementation of rnorm(10000, 10, 2) in R. It is presented as minimal demonstration; for a non-trivial application, see (LeBauer et al. 2013).

r.distn <- data.frame(distn = "norm", parama = 10, paramb = 2)
bugs.distn <- r2bugs.distributions(r.distn)
   
sample.bugs.distn <- function(prior = data.frame(distn = "norm", parama = 0, 
                                paramb = 1), n = 10000) {
  require(rjags)
  model.string <- paste0(
    "model{Y ~ d", prior$distn, 
    "(", prior$parama, 
    ## chisqr has only one parameter
    ifelse(prior$distn == "chisqr", "", paste0(", ", prior$paramb)), ");",
    "a <- x}"
  )   
  ## trick JAGS into running without data  
  writeLines(model.string, con = "test.bug")
  j.model  <- jags.model(file = "test.bug", data = list(x = 1))
  mcmc.object <- window(
    coda.samples(
      model = j.model, variable.names = c('Y'), 
      n.iter = n * 4, thin = 2),
    start = n)
  Y <- sample(as.matrix(mcmc.object)[,"Y"], n)
}
X <- sample.bugs.distn(bugs.distn)

4 Acknowlegements

This collaboration began on the Cross Validated statistical forum (http://stats.stackexchange.com/q/5543/1381). Funding was provided to DSL and MCD by the Energy Biosciences Institute.

CRAN packages used

rjags, R2WinBUGS

CRAN Task Views implied by cited packages

Bayesian, Cluster, GraphicalModels, MixedModels

Note

This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.