| find_lambda_one_d_rcpp {rust} | R Documentation |
Selecting the Box-Cox parameter in the 1D case using Rcpp
Description
Finds a value of the Box-Cox transformation parameter lambda for which
the (positive univariate) random variable with log-density
\log f has a density closer to that of a Gaussian random
variable. Works by estimating a set of quantiles of the distribution implied
by \log f and treating those quantiles as data in a standard
Box-Cox analysis. In the following we use theta (\theta) to
denote the argument of \log f on the original scale and
phi (\phi) on the Box-Cox transformed scale.
Usage
find_lambda_one_d_rcpp(
logf,
...,
ep_bc = 1e-04,
min_phi = ep_bc,
max_phi = 10,
num = 1001L,
xdiv = 100,
probs = seq(0.01, 0.99, by = 0.01),
lambda_range = c(-3, 3),
phi_to_theta = NULL,
log_j = NULL,
user_args = list()
)
Arguments
logf |
A pointer to a compiled C++ function returning the log
of the target density |
... |
further arguments to be passed to |
ep_bc |
A (positive) numeric scalar. Smallest possible value of
|
min_phi, max_phi |
Numeric scalars. Smallest and largest values
of |
num |
A numeric scalar. Number of values at which to evaluate
|
xdiv |
A numeric scalar. Only values of |
probs |
A numeric scalar. Probabilities at which to estimate the quantiles of that will be used as data to find lambda. |
lambda_range |
A numeric vector of length 2. Range of lambda over which to optimise. |
phi_to_theta |
A pointer to a compiled C++ function returning
(the inverse) of the transformation from |
log_j |
A pointer to a compiled C++ function returning the log of
the Jacobian of the transformation from |
user_args |
A list of numeric components providing arguments to
the user-supplied functions |
Details
The general idea is to estimate quantiles of f corresponding
to a set of equally-spaced probabilities in probs and to use these
estimated quantiles as data in a standard estimation of the Box-Cox
transformation parameter lambda.
The density f is first evaluated at num points equally spaced
over the interval (min_phi, max_phi). The continuous
density f is approximated by attaching trapezium-rule estimates of
probabilities to the midpoints of the intervals between the points. After
standardizing to account for the fact that f may not be normalized,
(min_phi, max_phi) is reset so that values with small
estimated probability (determined by xdiv) are excluded and the
procedure is repeated on this new range. Then the required quantiles are
estimated by inferring them from a weighted empirical distribution
function based on treating the midpoints as data and the estimated
probabilities at the midpoints as weights.
Value
A list containing the following components
lambda |
A numeric scalar. The value of |
gm |
A numeric scalar. Box-cox scaling parameter, estimated by the geometric mean of the quantiles used in the optimisation to find the value of lambda. |
init_psi |
A numeric scalar. An initial estimate of the mode of the Box-Cox transformed density |
sd_psi |
A numeric scalar. Estimates of the marginal standard deviations of the Box-Cox transformed variables. |
phi_to_theta |
as detailed above (only if |
log_j |
as detailed above (only if |
user_args |
as detailed above (only if |
References
Box, G. and Cox, D. R. (1964) An Analysis of Transformations. Journal of the Royal Statistical Society. Series B (Methodological), 26(2), 211-252.
Andrews, D. F. and Gnanadesikan, R. and Warner, J. L. (1971) Transformations of Multivariate Data, Biometrics, 27(4).
Eddelbuettel, D. and Francois, R. (2011). Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40(8), 1-18. doi:10.18637/jss.v040.i08
Eddelbuettel, D. (2013). Seamless R and C++ Integration with Rcpp, Springer, New York. ISBN 978-1-4614-6867-7.
See Also
ru_rcpp to perform ratio-of-uniforms sampling.
find_lambda_rcpp to produce (somewhat) automatically
a list for the argument lambda of ru for any value of
d.
Examples
# Log-normal density ===================
# Note: the default value of max_phi = 10 is OK here but this will not
# always be the case.
ptr_lnorm <- create_xptr("logdlnorm")
mu <- 0
sigma <- 1
lambda <- find_lambda_one_d_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma)
lambda
x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, log = TRUE, d = 1,
n = 1000, trans = "BC", lambda = lambda)
# Gamma density ===================
alpha <- 1
# Choose a sensible value of max_phi
max_phi <- qgamma(0.999, shape = alpha)
# [I appreciate that typically the quantile function won't be available.
# In practice the value of lambda chosen is quite insensitive to the choice
# of max_phi, provided that max_phi is not far too large or far too small.]
ptr_gam <- create_xptr("logdgamma")
lambda <- find_lambda_one_d_rcpp(logf = ptr_gam, alpha = alpha,
max_phi = max_phi)
lambda
x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = 1000, trans = "BC",
lambda = lambda)
alpha <- 0.1
# NB. for alpha < 1 the gamma(alpha, beta) density is not bounded
# So the ratio-of-uniforms emthod can't be used but it may work after a
# Box-Cox transformation.
# find_lambda_one_d() works much better than find_lambda() here.
max_phi <- qgamma(0.999, shape = alpha)
lambda <- find_lambda_one_d_rcpp(logf = ptr_gam, alpha = alpha,
max_phi = max_phi)
lambda
x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = 1000, trans = "BC",
lambda = lambda)
plot(x)
plot(x, ru_scale = TRUE)