nimbleNoBounds (version 1.0.2)

dLogitBeta: Logit transformed beta distribution.

Description

Logit transformation of beta random variable to the real line.

Usage

dLogitBeta(x, shape1 = 1, shape2 = 1, log = 0)
rLogitBeta(n = 1, shape1 = 1, shape2 = 1)

Value

density or log-density of beta distributions transformed to real line via logit function.

Arguments

x

A continuous random variable on the real line, where y=ilogit(x) and y ~ dbeta(shape1, shape2).

shape1

non-negative parameter of the Beta distribution.

shape2

non-negative parameter of the Beta distribution.

log

logical flag. Returns log-density if TRUE.

n

Number of random variables. Currently limited to 1, as is common in nimble. See ?replicate for an alternative.

Author

David R.J. Pleydell

Examples

## Create a beta random variable, and transform it to the logit scale
n      = 100000
sh1    = 1
sh2    = 11
y      = rbeta(n=n, sh1, sh2)
x      = logit(y)

## Plot histograms of the two random variables
oldpar <- par()
par(mfrow=n2mfrow(2))
## Plot 1
hist(x, n=100, freq=FALSE)
curve(dLogitBeta(x, sh1, sh2), -15, 4, n=1001, col="red", add=TRUE, lwd=3)
## Plot 2: back-transformed
xNew = replicate(n=n, rLogitBeta(n=1, sh1, sh2))
yNew   = ilogit(xNew)
hist(yNew, n=100, freq=FALSE, xlab="exp(x)")
curve(dbeta(x, sh1, sh2), 0, 1, n=1001, col="red", lwd=3, add=TRUE)
par(oldpar)

## Create a NIMBLE model that uses this transformed distribution
code = nimbleCode({
  x ~ dLogitBeta(sh1, sh2)
})

# \donttest{
## Build & compile the model
const = list(sh1=sh1, sh2=sh2)
modelR = nimbleModel(code=code, const=const)
modelC = compileNimble(modelR)

## Configure, build and compile an MCMC
conf  = configureMCMC(modelC)
mcmc  = buildMCMC(conf=conf)
cMcmc = compileNimble(mcmc)

## Run the MCMC
x = as.vector(runMCMC(mcmc=cMcmc, niter=50000))

## Plot MCMC output
oldpar <- par()
par(mfrow=n2mfrow(3))
## Plot 1: MCMC trajectory
plot(x, typ="l")
## Plot 2: taget density on unbounded sampling scale
hist(x, n=100, freq=FALSE, xlab="x = logit(y)")
curve(dLogitBeta(x, sh1, sh2), -15, 5, n=1001, col="red", lwd=3, add=TRUE)
## Plot 3: taget density on bounded scale
hist(ilogit(x), n=100, freq=FALSE, xlab="y")
curve(dbeta(x, sh1, sh2), 0, 1, n=1001, col="red", lwd=3, add=TRUE)
par(oldpar)
# }