dLogChisq {nimbleNoBounds}

Log transformed chi-squared distribution.

Description

dLogChisq and rLogChisq provide a log-transformed chi-squared distribution.

Usage

dLogChisq(x, df = 1, log = 0)

rLogChisq(n = 1, df = 1)

Arguments

x

A continuous random variable on the real line, where y=exp(x) and y ~ dchisq(df).

df

df parameter of y ~ dchisq(df).

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.

Value

The density or log density of x, such that x = log(y) and y ~ dchisq(df).

Author(s)

David R.J. Pleydell

Examples

## Create a Chi-squared random variable, and transform it to the log scale
n  = 100000
df = 3
y  = rchisq(n=n, df=df)
x  = log(y)

## Plot histograms of the two random variables
oldpar <- par()
par(mfrow=n2mfrow(2))
## Plot 1
hist(x, n=100, freq=FALSE, main="", xlab= "x = log(y)")
curve(dLogChisq(x, df=df), -7, 4, n=1001, col="red", add=TRUE, lwd=3)
## Plot 2: back-transformed
xNew = replicate(n=n, rLogChisq(n=1, df=df))
yNew = exp(xNew)
hist(yNew, n=100, freq=FALSE, xlab="y = exp(x)", main="")
curve(dchisq(x, df=df), 0, 30, n=1001, col="red", lwd=3, add=TRUE)
par(oldpar)

## Create a NIMBLE model that uses this transformed distribution
code = nimbleCode({
  x ~ dLogChisq(df = 0.5)
  y <- exp(x)
})


## Build & compile the model
modelR = nimbleModel(code=code)
modelC = compileNimble(modelR)

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

## Run the MCMC
samps = runMCMC(mcmc=cMcmc, niter=50000)
x     = samps[,"x"]
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, main="Histogram of MCMC samples", xlab="x = log(y)")
curve(dLogChisq(x, df=0.5), -55, 5, 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="Back-transformed MCMC samples")
curve(dchisq(x, df=0.5), 0, 20, n=1001, col="red", lwd=3, add=TRUE)
par(oldpar)