morph(b, r, p, center)
morph.identity()p is
specified, r defaults to 0.r is specified, p defaults to 3.center is a vector it should be the same length of the state
of the Markov chain, center defaults to 0outfun(f), a function that operates on functions.outfun(f)returns the functionfunction(state, ...)
f(inverse(state), ...).inverse, the inverse transformation function.transform, the transformation function.lud, a function that operates on functions. As input,ludtakes 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 thetransformfunction.lud(f) = function(state, ...)
f(inverse(state), ...) + log.jacobian(state), wherelog.jacobianrepresents the function that calculate the log
Jacobian of the transformation.log.jacobianis not 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.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$):
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:
randpis specified, andbis not specified, then$$f^{-1}(X) = g_1(|X|)
\frac{X}{|X|}.$$If onlyris specified,p = 3is used. If onlypis specified,r = 0is used.bis specified, then$$f^{-1}(X) = g_2(|X|)
\frac{X}{|X|}.$$randpis specified, andbis
also specified, then$$f^{-1}(X) = g_2(g_1(|X|))
\frac{X}{|X|}.$$morph.metrop# 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)))
curve(d.transformed, from=-3, to=3, ylab="Induced Density")Run the code above in your browser using DataLab