Last data update: 2014.03.03

R: Variable Transformation
morphR 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(|X|) * X / |X|

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

f^{-1}(Y) = g^{-1}(|Y|) * Y / |Y|.

The parameters r and p are used to define the function

g1(x) = x + (x-r)^p * I(x > r)

where I(•) 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

g2(x) = (exp(b * x) - exp(1) / 3) * I(x > 1 / b) + (x^3 * b^3 exp(1) / 6 + x * b * exp(1) / 2) * I(x <= 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) = g1(|X|) * 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) = g2(|X|) * X / |X|.

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

    f^{-1}(X) = g2(g1(|X|)) * 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.

See Also

morph.metrop

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)

Results