morph: Variable Transformation
In mcmc: Markov Chain Monte Carlo

morph

R Documentation

Variable Transformation

Description

Utility functions for variable transformation.

Usage

morph(b, r, p, center)
morph.identity()

Arguments

`b`	Positive real number. May be missing.
`r`	Non-negative real number. May be missing. If `p` is specified, `r` defaults to 0.
`p`	Real number strictly greater than 2. May be missing. If `r` is specified, `p` defaults to 3.
`center`	Real scalar or vector. May be missing. If `center` is a vector it should be the same length of the state of the Markov chain, `center` defaults to 0

Details

The morph function facilitates using variable transformations by providing functions to (using X for the original random variable with the pdf f_X, and Y for the transformed random variable with the pdf f_Y):

Calculate the log unnormalized probability density for Y induced by the transformation.
Transform an arbitrary function of X to a function of Y.
Transform values of X to values of Y.
Transform values of Y to values of X (the inverse transformation).

for a select few transformations.

morph.identity implements the identity transformation, Y=X.

The parameters r, p, b and center specify the transformation function. In all cases, center gives the center of the transformation, which is the value c in the equation

Y = f(X - c).

If no parameters are specified, the identity transformation, Y=X, is used.

The parameters r, p and b specify a function g, which is a monotonically increasing bijection from the non-negative reals to the non-negative reals. Then

f(X) = g\bigl(|X|\bigr) \frac{X}{|X|}

where |X| represents the Euclidean norm of the vector X. The inverse function is given by

f^{-1}(Y) = g^{-1}\bigl(|Y|\bigr) \frac{Y}{|Y|}.

The parameters r and p are used to define the function

g_1(x) = x + (x-r)^p I(x > r)

where I( \cdot ) is the indicator function. We require that r is non-negative and p is strictly greater than 2. The parameter b is used to define the function

g_2(x) = \bigl(e^{bx} - e / 3\bigr) I(x > \frac{1}{b}) + \bigl(x^3 b^3 e / 6 + x b e / 2\bigr) I(x \leq \frac{1}{b})

We require that b is positive.

The parameters r, p and b specify f^{-1} in the following manner:

If one or both of r and p is specified, and b is not specified, then

f^{-1}(X) = g_1(|X|) \frac{X}{|X|}.

If only r is specified, p = 3 is used. If only p is specified, r = 0 is used.
If only b is specified, then

f^{-1}(X) = g_2(|X|) \frac{X}{|X|}.
If one or both of r and p is specified, and b is also specified, then

f^{-1}(X) = g_2(g_1(|X|)) \frac{X}{|X|}.

Value

a list containing the functions

outfun(f), a function that operates on functions. outfun(f) returns the function function(state, ...) f(inverse(state), ...).
inverse, the inverse transformation function.
transform, the transformation function.
lud, a function that operates on functions. As input, lud takes a function that calculates a log unnormalized probability density, and returns a function that calculates the log unnormalized density by transforming a random variable using the transform function. lud(f) = function(state, ...) f(inverse(state), ...) + log.jacobian(state), where log.jacobian represents the function that calculate the log Jacobian of the transformation. log.jacobian is not returned.

Warning

The equations for the returned transform function (see below) do not have a general analytical solution when p is not equal to 3. This implementation uses numerical approximation to calculate transform when p is not equal to 3. If computation speed is a factor, it is advisable to use p=3. This is not a factor when using morph.metrop, as transform is only called once during setup, and not at all while running the Markov chain.

Examples


# use an exponential transformation, centered at 100.
b1 <- morph(b=1, center=100)
# original log unnormalized density is from a t distribution with 3
# degrees of freedom, centered at 100.
lud.transformed <- b1$lud(function(x) dt(x - 100, df=3, log=TRUE))
d.transformed <- Vectorize(function(x) exp(lud.transformed(x)))
## Not run: 
curve(d.transformed, from=-3, to=3, ylab="Induced Density")

## End(Not run)

mcmc documentation built on Nov. 17, 2023, 1:06 a.m.