nimbleNoBounds (version 1.0.2)

dLogWeib: Log transformed Weibull distribution.

Description

dLogWeib and rLogWeib provide a log-transformed Weibull distribution.

Usage

dLogWeib(x, shape = 1, scale = 1, log = 0)
rLogWeib(n = 1, shape = 1, scale = 1)

Value

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

Arguments

x

A continuous random variable on the real line, where y=exp(x) and y ~ dweib(shape,scale).

shape

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

scale

Scale parameter of y ~ dweib(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 a Weibull random variable, and transform it to the log scale
n      = 100000
shape = 2
scale = 2
y      = rweibull(n=n, shape=shape, scale=scale)
x      = log(y)

## Plot histograms of the two random variables
oldpar <- par()
par(mfrow=n2mfrow(2))
## Plot 1
hist(x, n=100, freq=FALSE, xlab="x = log(y)",
     main="Histogram of log-transformed random numbers from rweibull.")
curve(dLogWeib(x, shape=shape, scale=scale), -6, 3, n=1001, col="red", add=TRUE, lwd=3)
## Plot 2: back-transformed
xNew = replicate(n=n, rLogWeib(n=1, shape=shape, scale=scale))
yNew = exp(xNew)
hist(yNew, n=100, freq=FALSE, xlab="y = exp(x)", main="Histogram of random numbers from rLogWeib.")
curve(dweibull(x, shape=shape, scale=scale), 0, 100, n=1001, col="red", lwd=3, add=TRUE)
par(oldpar)

## Create a NIMBLE model that uses this transformed distribution
code = nimbleCode({
  log(y) ~ dLogWeib(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)
mcmc  = buildMCMC(conf=conf)
cMcmc = compileNimble(mcmc)

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

## 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, main="Histogram of MCMC samples", xlab="x = log(y)")
curve(dLogWeib(x, shape=shape, scale=scale), -5, 3, n=1001, col="red", lwd=3, add=TRUE)
## Plot 3: taget density on bounded scale
hist(y, n=100, freq=FALSE, xlab="y = exp(x)",
     main="Histogram of back-transformed MCMC samples")
curve(dweibull(x, shape=shape, scale=scale), 0, 8, n=1001, col="red", lwd=3, add=TRUE)
par(oldpar)
# }