ru: Generalized ratio-of-uniforms sampling

Description Usage Arguments Details Value References See Also Examples

View source: R/ru_sampling.R

Description

Uses the generalized ratio-of-uniforms method to simulate from a distribution with log-density log f (up to an additive constant). The density f must be bounded, perhaps after a transformation of variable.

Usage

1
2
3
4
5
6
7
8
9
ru(logf, ..., n = 1, d = 1, init = NULL, trans = c("none", "BC",
  "user"), phi_to_theta = NULL, log_j = NULL, user_args = list(),
  lambda = rep(1L, d), lambda_tol = 1e-06, gm = NULL, rotate = ifelse(d
  == 1, FALSE, TRUE), lower = rep(-Inf, d), upper = rep(Inf, d), r = 1/2,
  ep = 0L, a_algor = if (d == 1) "nlminb" else "optim",
  b_algor = c("nlminb", "optim"), a_method = c("Nelder-Mead", "BFGS", "CG",
  "L-BFGS-B", "SANN", "Brent"), b_method = c("Nelder-Mead", "BFGS", "CG",
  "L-BFGS-B", "SANN", "Brent"), a_control = list(), b_control = list(),
  var_names = NULL)

Arguments

logf

A function returning the log of the target density f. This function should return -Inf when the density is zero.

...

Further arguments to be passed to logf and related functions.

n

A numeric scalar. Number of simulated values required.

d

A numeric scalar. Dimension of f.

init

A numeric vector. Initial estimates of the mode of logf. If trans="BC" or trans = "user" this is after Box-Cox transformation or user-defined transformation, but before any rotation of axes.

trans

A character scalar. "none" for no transformation, "BC" for Box-Cox transformation, "user" for a user-defined transformation. If trans = "user" then the transformation should be specified using phi_to_theta and log_j and user_args may be used to pass arguments to phi_to_theta and log_j.

phi_to_theta

A function returning (the inverse) of the transformation from theta to phi used to ensure positivity of phi prior to Box-Cox transformation. The argument is phi and the returned value is theta. If phi_to_theta is undefined at the input value then the function should return NA.

log_j

A function returning the log of the Jacobian of the transformation from theta to phi, i.e. based on derivatives of phi with respect to theta. Takes theta as its argument.

user_args

A list of numeric components. If trans = ``user'' then user_args is a list providing arguments to the user-supplied functions phi_to_theta and log_j.

lambda

Either

  • A numeric vector. Box-Cox transformation parameters, or

  • A list with components

    lambda

    A numeric vector. Box-Cox parameters (required).

    gm

    A numeric vector. Box-cox scaling parameters (optional). If supplied this overrides any gm supplied by the individual gm argument described below.

    init_psi

    A numeric vector. Initial estimate of mode after Box-Cox transformation (optional).

    sd_psi

    A numeric vector. Estimates of the marginal standard deviations of the Box-Cox transformed variables (optional).

    phi_to_theta

    as above (optional).

    log_j

    as above (optional).

    This list may be created using find_lambda_one_d (for d = 1) or find_lambda (for any d).

lambda_tol

A numeric scalar. Any values in lambda that are less than lambda_tol in magnitude are set to zero.

gm

A numeric vector. Box-cox scaling parameters (optional). If lambda$gm is supplied in input list lambda then lambda$gm is used, not gm.

rotate

A logical scalar. If TRUE (d > 1 only) use Choleski rotation. If d = 1 and rotate = TRUE then rotate will be set to FALSE with a warning.

lower, upper

Numeric vectors. Lower/upper bounds on the arguments of the function after any transformation from theta to phi implied by the inverse of phi_to_theta. If rotate = FALSE these are used in all of the optimizations used to construct the bounding box. If rotate = TRUE then they are use only in the first optimisation to maximise the target density.' If trans = "BC" components of lower that are negative are set to zero without warning and the bounds implied after the Box-Cox transformation are calculated inside ru.

r

A numeric scalar. Parameter of generalized ratio-of-uniforms.

ep

A numeric scalar. Controls initial estimates for optimizations to find the b-bounding box parameters. The default (ep=0) corresponds to starting at the mode of logf small positive values of ep move the constrained variable slightly away from the mode in the correct direction. If ep is negative its absolute value is used, with no warning given.

a_algor, b_algor

Character scalars. Either "nlminb" or "optim". Respective optimization algorithms used to find a(r) and (bi-(r), bi+(r)).

a_method, b_method

Character scalars. Respective methods used by optim to find a(r) and (bi-(r), bi+(r)). Only used if optim is the chosen algorithm. If d = 1 then a_method and b_method are set to "Brent" without warning.

a_control, b_control

Lists of control arguments to optim or nlminb to find a(r) and (bi-(r), bi+(r)) respectively.

var_names

A character vector. Names to give to the column(s) of the simulated values.

Details

If trans = "none" and rotate = FALSE then ru implements the (multivariate) generalized ratio of uniforms method described in Wakefield, Gelfand and Smith (1991) using a target density whose mode is relocated to the origin (‘mode relocation’) in the hope of increasing efficiency.

If trans = "BC" then marginal Box-Cox transformations of each of the d variables is performed, with parameters supplied in lambda. The function phi_to_theta may be used, if necessary, to ensure positivity of the variables prior to Box-Cox transformation.

If trans = "user" then the function phi_to_theta enables the user to specify their own transformation.

In all cases the mode of the target function is relocated to the origin after any user-supplied transformation and/or Box-Cox transformation.

If d is greater than one and rotate = TRUE then a rotation of the variable axes is performed after mode relocation. The rotation is based on the Choleski decomposition (see chol) of the estimated Hessian (computed using optimHess of the negated log-density after any user-supplied transformation or Box-Cox transformation. If any of the eigenvalues of the estimated Hessian are non-positive (which may indicate that the estimated mode of logf is close to a variable boundary) then rotate is set to FALSE with a warning. A warning is also given if this happens when d = 1.

The default value of the tuning parameter r is 1/2, which is likely to be close to optimal in many cases, particularly if trans = "BC".

See vignette("rust-vignette", package = "rust") for full details.

Value

An object of class "ru" is a list containing the following components:

sim_vals

An n by d matrix of simulated values.

box

A (2 * d + 1) by d + 2 matrix of ratio-of-uniforms bounding box information, with row names indicating the box parameter. The columns contain

column 1

values of box parameters.

columns 2 to (2+d-1)

values of variables at which these box parameters are obtained.

column 2+d

convergence indicators.

Scaling of f within ru and relocation of the mode to the origin means that the first row of box will always be c(1, rep(0, d)).

pa

A numeric scalar. An estimate of the probability of acceptance.

d

A numeric scalar. The dimension of logf.

logf

A function. logf supplied by the user, but with f scaled by the maximum of the target density used in the ratio-of-uniforms method (i.e. logf_rho), to avoid numerical problems in contouring fin plot.ru when d = 2.

logf_rho

A function. The target function actually used in the ratio-of-uniforms algorithm.

sim_vals_rho

An n by d matrix of values simulated from the function used in the ratio-of-uniforms algorithm.

logf_args

A list of further arguments to logf.

f_mode

The estimated mode of the target density f, after any Box-Cox transformation and/or user supplied transformation, but before mode relocation.

References

Wakefield, J. C., Gelfand, A. E. and Smith, A. F. M. (1991) Efficient generation of random variates via the ratio-of-uniforms method. Statistics and Computing (1991), 1, 129-133. http://dx.doi.org/10.1007/BF01889987.

See Also

ru_rcpp for a version of ru that uses the Rcpp package to improve efficiency.

summary.ru for summaries of the simulated values and properties of the ratio-of-uniforms algorithm.

plot.ru for a diagnostic plot (for d = 1 and d = 2 only).

find_lambda_one_d to produce (somewhat) automatically a list for the argument lambda of ru for the d = 1 case.

find_lambda to produce (somewhat) automatically a list for the argument lambda of ru for any value of d.

optim for choices of the arguments a_method, b_method, a_control and b_control.

nlminb for choices of the arguments a_control and b_control.

optimHess for Hessian estimation.

chol for the Choleski decomposition.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
# Normal density ===================

# One-dimensional standard normal ----------------
x <- ru(logf = function(x) -x ^ 2 / 2, d = 1, n = 1000, init = 0.1)

# Two-dimensional standard normal ----------------
x <- ru(logf = function(x) -(x[1]^2 + x[2]^2) / 2, d = 2, n = 1000,
        init = c(0, 0))

# Two-dimensional normal with positive association ----------------
rho <- 0.9
covmat <- matrix(c(1, rho, rho, 1), 2, 2)
log_dmvnorm <- function(x, mean = rep(0, d), sigma = diag(d)) {
  x <- matrix(x, ncol = length(x))
  d <- ncol(x)
  - 0.5 * (x - mean) %*% solve(sigma) %*% t(x - mean)
}

# No rotation.
x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0),
        rotate = FALSE)

# With rotation.
x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0))

# three-dimensional normal with positive association ----------------
covmat <- matrix(rho, 3, 3) + diag(1 - rho, 3)

# No rotation.  Slow !
x <- ru(logf = log_dmvnorm, sigma = covmat, d = 3, n = 1000,
        init = c(0, 0, 0), rotate = FALSE)

# With rotation.
x <- ru(logf = log_dmvnorm, sigma = covmat, d = 3, n = 1000,
        init = c(0, 0, 0))

# Log-normal density ===================

# Sampling on original scale ----------------
x <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, lower = 0, init = 1)

# Box-Cox transform with lambda = 0 ----------------
lambda <- 0
x <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, lower = 0, init = 0.1,
        trans = "BC", lambda = lambda)

# Equivalently, we could use trans = "user" and supply the (inverse) Box-Cox
# transformation and the log-Jacobian by hand
x <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, init = 0.1,
        trans = "user", phi_to_theta = function(x) exp(x),
        log_j = function(x) -log(x))

# Gamma(alpha, 1) density ===================

# Note: the gamma density in unbounded when its shape parameter is < 1.
# Therefore, we can only use trans="none" if the shape parameter is >= 1.

# Sampling on original scale ----------------

alpha <- 10
x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000,
        lower = 0, init = alpha)

alpha <- 1
x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000,
        lower = 0, init = alpha)

# Box-Cox transform with lambda = 1/3 works well for shape >= 1. -----------

alpha <- 1
x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000,
        trans = "BC", lambda = 1/3, init = alpha)
summary(x)

# Equivalently, we could use trans = "user" and supply the (inverse) Box-Cox
# transformation and the log-Jacobian by hand

# Note: when phi_to_theta is undefined at x this function returns NA
phi_to_theta  <- function(x, lambda) {
  ifelse(x * lambda + 1 > 0, (x * lambda + 1) ^ (1 / lambda), NA)
}
log_j <- function(x, lambda) (lambda - 1) * log(x)
lambda <- 1/3
x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000,
        trans = "user", phi_to_theta = phi_to_theta, log_j = log_j,
        user_args = list(lambda = lambda), init = alpha)
summary(x)

## Not run: 
# Generalized Pareto posterior distribution ===================

# Sample data from a GP(sigma, xi) distribution
gpd_data <- rgpd(m = 100, xi = -0.5, sigma = 1)
# Calculate summary statistics for use in the log-likelihood
ss <- gpd_sum_stats(gpd_data)
# Calculate an initial estimate
init <- c(mean(gpd_data), 0)

# Mode relocation only ----------------
n <- 1000
x1 <- ru(logf = gpd_logpost, ss = ss, d = 2, n = n, init = init,
         lower = c(0, -Inf), rotate = FALSE)
plot(x1, xlab = "sigma", ylab = "xi")
# Parameter constraint line xi > -sigma/max(data)
# [This may not appear if the sample is far from the constraint.]
abline(a = 0, b = -1 / ss$xm)
summary(x1)

# Rotation of axes plus mode relocation ----------------
x2 <- ru(logf = gpd_logpost, ss = ss, d = 2, n = n, init = init,
         lower = c(0, -Inf))
plot(x2, xlab = "sigma", ylab = "xi")
abline(a = 0, b = -1 / ss$xm)
summary(x2)

# Cauchy ========================

# The bounding box cannot be constructed if r < 1.  For r = 1 the
# bounding box parameters b1-(r) and b1+(r) are attained in the limits
# as x decreases/increases to infinity respectively.  This is fine in
# theory but using r > 1 avoids this problem and the largest probability
# of acceptance is obtained for r approximately equal to 1.26.

res <- ru(logf = dcauchy, log = TRUE, init = 0, r = 1.26, n = 1000)

# Half-Cauchy ===================

log_halfcauchy <- function(x) {
  return(ifelse(x < 0, -Inf, dcauchy(x, log = TRUE)))
}

# Like the Cauchy case the bounding box cannot be constructed if r < 1.
# We could use r > 1 but the mode is on the edge of the support of the
# density so as an alternative we use a log transformation.

x <- ru(logf = log_halfcauchy, init = 0, trans = "BC", lambda = 0, n = 1000)
x$pa
plot(x, ru_scale = TRUE)

# Example 4 from Wakefield et al. (1991) ===================

# Bivariate normal x bivariate student-t
log_norm_t <- function(x, mean = rep(0, d), sigma1 = diag(d), sigma2 = diag(d)) {
  x <- matrix(x, ncol = length(x))
  d <- ncol(x)
  log_h1 <- -0.5 * (x - mean) %*% solve(sigma1) %*% t(x - mean)
  log_h2 <- -2 * log(1 + 0.5 * x %*% solve(sigma2) %*% t(x))
  return(log_h1 + log_h2)
}

rho <- 0.9
covmat <- matrix(c(1, rho, rho, 1), 2, 2)
y <- c(0, 0)

# Case in the top right corner of Table 3
x <- ru(logf = log_norm_t, mean = y, sigma1 = covmat, sigma2 = covmat,
  d = 2, n = 10000, init = y, rotate = FALSE)
x$pa

# Rotation increases the probability of acceptance
x <- ru(logf = log_norm_t, mean = y, sigma1 = covmat, sigma2 = covmat,
  d = 2, n = 10000, init = y, rotate = TRUE)
x$pa

## End(Not run)

rust documentation built on May 16, 2018, 5:06 p.m.