nimbleNoBounds (version 1.0.2)

dLogInvgamma: Log transformed inverse-gamma distribution.

Description

dLogInvgamma and rLogInvgamma provide a log-transformed inverse gamma distribution.

Usage

dLogInvgamma(x, shape, scale = 1, log = 0)
rLogInvgamma(n = 1, shape, scale = 1)

Value

The density or log density of x, such that x = log(y) and y ~ dinvgamma(shape,scale).

Arguments

x

A continuous random variable on the real line. Let y=exp(x). Then y ~ dinvgamma(shape,scale).

shape

Shape parameter of y ~ dinvgamma(shape,scale).

scale

Scale parameter of y ~ dinvgamma(shape,scale).

log

Logical flag to toggle returning the log density.

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 an inverse gamma random variable, and transform it to the log scale
n      = 100000
shape = 2.5
scale = 0.01
y      = rinvgamma(n=n, shape=shape, scale=scale)
x      = log(y)

## Plot histograms of the two random variables
oldpar <- par()
par(mfrow=n2mfrow(4))
## Plot 1
hist(y, n=100, freq=FALSE, xlab="y")
curve(dinvgamma(x, shape=shape, scale=scale), 0, 1.0, n=5001, col="red", add=TRUE, lwd=2)
## Plot 2
hist(x, n=100, freq=FALSE)
curve(dLogInvgamma(x, shape=shape, scale=scale), -8, 1, n=1001, col="red", add=TRUE, lwd=3)
## Plot 3: back-transformed
z = rgamma(n=n, shape=shape, scale=1/scale)
hist(1/z, n=100, freq=FALSE, xlab="y")
curve(dinvgamma(x, shape=shape, scale=scale), 0, 1, n=5001, col="red", lwd=3, add=TRUE)
## Plot 4: back-transformed
xNew = replicate(n=n, rLogInvgamma(n=1, shape=shape, scale=scale))
yNew = exp(xNew)
hist(yNew, n=100, freq=FALSE, xlab="exp(x)")
curve(dinvgamma(x, shape=shape, scale=scale), 0, 1, n=5001, col="red", lwd=3, add=TRUE)
par(oldpar)

## Create a NIMBLE model that uses this transformed distribution
code = nimbleCode({
  log(y) ~ dLogInvgamma(shape=shape, scale=scale)
})

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

## Configure, build and compile an MCMC
conf  = configureMCMC(modelC, monitors=c("log_y", "y"))
mcmc  = buildMCMC(conf=conf)
cMcmc = compileNimble(mcmc)

## Run the MCMC
samps = runMCMC(mcmc=cMcmc, niter=50000)
x = samps[,"log_y"]
y = samps[,"y"]


## 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)
curve(dLogInvgamma(x, shape=shape, scale=scale), -10, 1, n=1001, col="red", lwd=3, add=TRUE)
## Plot 3: taget density on bounded scale
curve(dinvgamma(x, shape=shape, scale=scale), xlab="y = exp(x)", 0, 0.5, n=1001, col="red", lwd=3)
hist(y, n=100, freq=FALSE, add=TRUE)
curve(dinvgamma(x, shape=shape, scale=scale), 0, 0.5, n=1001, col="red", add=TRUE, lwd=3)
par(oldpar)
# }