Movatterモバイル変換

rust: Ratio-of-Uniforms Simulation with Transformation

Description

Uses the multivariate generalized ratio-of-uniforms method to simulate from adistribution with log-densitylogf (up to an additive constant).logf must be bounded, perhaps after a transformation of variable.

Details

The main functions in the rust package areru andru_rcpp, which implement the generalized ratio-of-uniformsalgorithm. The latter uses the Rcpp package to improve efficiency.Also provided are two functions,find_lambda andfind_lambda_one_d, that may beused to set a suitable value for the parameterlambda if Box-Coxtransformation is used prior to simulation.Ifru_rcpp is used the equivalent functions arefind_lambda_rcpp andfind_lambda_one_d_rcppBasicplot andsummary methods are also provided.

See the following package vignettes for information:

• Introducing rust orvignette("rust-a-vignette", package = "rust").

• When can rust be used? orvignette("rust-b-when-to-use-vignette", package = "rust").

• Rusting faster: Simulation using Rcpp orvignette("rust-c-using-rcpp-vignette", package = "rust").

Author(s)

Maintainer: Paul J. Northropp.northrop@ucl.ac.uk [copyright holder]

References

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

Box, G. and Cox, D. R. (1964) An Analysis of Transformations.Journal of the Royal Statistical Society. Series B (Methodological), 26(2),211-252.

Eddelbuettel, D. and Francois, R. (2011). Rcpp: Seamless R andC++ Integration. Journal of Statistical Software, 40(8), 1-18.doi:10.18637/jss.v040.i08.

Eddelbuettel, D. (2013) Seamless R and C++ Integration withRcpp. Springer, New York. ISBN 978-1-4614-6867-7.

See Also

ru andru_rcpp to performratio-of-uniforms sampling.

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

plot.ru for a diagnostic plot.

find_lambda_one_d andfind_lambda_one_d_rcpp to produce (somewhat) automaticallya list for the argumentlambda ofru for thed = 1 case.

find_lambda andfind_lambda_rcppto produce (somewhat) automaticallya list for the argumentlambda ofru for any value ofd.

Create external pointer to a C++ function forlog_j

Description

Create external pointer to a C++ function forlog_j

Usage

create_log_j_xptr(fstr)

Arguments

Details

See theRusting faster: Simulation using Rcpp vignette.

Examples

See the examples inru_rcpp.

Create external pointer to a C++ function forphi_to_theta

Description

Create external pointer to a C++ function forphi_to_theta

Usage

create_phi_to_theta_xptr(fstr)

Arguments

Details

See theRusting faster: Simulation using Rcpp vignette.

Examples

See the examples inru_rcpp.

Create external pointer to a C++ function forlogf

Description

Create external pointer to a C++ function forlogf

Usage

create_xptr(fstr)

Arguments

Details

See theRusting faster: Simulation using Rcpp vignette.

Examples

See the examples inru_rcpp.

Selecting the Box-Cox parameter for general d

Description

Finds a value of the Box-Cox transformation parameter lambda for whichthe (positive) random variable with log-density\log f has adensity closer to that of a Gaussian random variable.In the following we usetheta (\theta) to denote the argumentof\log f on the original scale andphi (\phi) onthe Box-Cox transformed scale.

Usage

find_lambda(  logf,  ...,  d = 1,  n_grid = NULL,  ep_bc = 1e-04,  min_phi = rep(ep_bc, d),  max_phi = rep(10, d),  which_lam = 1:d,  lambda_range = c(-3, 3),  init_lambda = NULL,  phi_to_theta = NULL,  log_j = NULL)

Arguments

Details

The general idea is to evaluate the densityf on ad-dimensional grid, withn_grid ordinates for each of thed variables.We treat each combination of the variables in the grid as a data pointand perform an estimation of the Box-Cox transformation parameterlambda, in which each data point is weighted by the densityat that point. The vectorsmin_phi andmax_phi define thelimits of the grid andwhich_lam can be used to specify that onlycertain components ofphi are to be transformed.

Value

A list containing the following components

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).

See Also

ru andru_rcpp to performratio-of-uniforms sampling.

find_lambda_one_d andfind_lambda_one_d_rcpp to produce (somewhat) automaticallya list for the argument lambda ofru/ru_rcpp for thed = 1 case.

find_lambda_rcpp for a version offind_lambda that uses the Rcpp package to improveefficiency.

Examples

# Log-normal density ===================# Note: the default value max_phi = 10 is OK here but this will not always# be the caselambda <- find_lambda(logf = dlnorm, log = TRUE)lambdax <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, trans = "BC",        lambda = lambda)# Gamma density ===================alpha <- 1#  Choose a sensible value of max_phimax_phi <- qgamma(0.999, shape = alpha)# [Of course, typically the quantile function won't be available.  However,# 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.]lambda <- find_lambda(logf = dgamma, shape = alpha, log = TRUE,                      max_phi = max_phi)lambdax <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000,        trans = "BC", lambda = lambda)# Generalized Pareto posterior distribution ===================# Sample data from a GP(sigma, xi) distributiongpd_data <- rgpd(m = 100, xi = -0.5, sigma = 1)# Calculate summary statistics for use in the log-likelihoodss <- gpd_sum_stats(gpd_data)# Calculate an initial estimateinit <- c(mean(gpd_data), 0)n <- 1000# Sample on original scale, with no rotation ----------------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)# Sample on original scale, with rotation ----------------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)# Sample on Box-Cox transformed scale ----------------# Find initial estimates for phi = (phi1, phi2),# where phi1 = sigma#   and phi2 = xi + sigma / max(x),# and ranges of phi1 and phi2 over over which to evaluate# the posterior to find a suitable value of lambda.temp <- do.call(gpd_init, ss)min_phi <- pmax(0, temp$init_phi - 2 * temp$se_phi)max_phi <- pmax(0, temp$init_phi + 2 * temp$se_phi)# Set phi_to_theta() that ensures positivity of phi# We use phi1 = sigma and phi2 = xi + sigma / max(data)phi_to_theta <- function(phi) c(phi[1], phi[2] - phi[1] / ss$xm)log_j <- function(x) 0lambda <- find_lambda(logf = gpd_logpost, ss = ss, d = 2, min_phi = min_phi,  max_phi = max_phi, phi_to_theta = phi_to_theta, log_j = log_j)lambda# Sample on Box-Cox transformed, without rotationx3 <- ru(logf = gpd_logpost, ss = ss, d = 2, n = n, trans = "BC",  lambda = lambda, rotate = FALSE)plot(x3, xlab = "sigma", ylab = "xi")abline(a = 0, b = -1 / ss$xm)summary(x3)# Sample on Box-Cox transformed, with rotationx4 <- ru(logf = gpd_logpost, ss = ss, d = 2, n = n, trans = "BC",  lambda = lambda)plot(x4, xlab = "sigma", ylab = "xi")abline(a = 0, b = -1 / ss$xm)summary(x4)def_par <- graphics::par(no.readonly = TRUE)par(mfrow = c(2,2), mar = c(4, 4, 1.5, 1))plot(x1, xlab = "sigma", ylab = "xi", ru_scale = TRUE,  main = "mode relocation")plot(x2, xlab = "sigma", ylab = "xi", ru_scale = TRUE,  main = "mode relocation and rotation")plot(x3, xlab = "sigma", ylab = "xi", ru_scale = TRUE,  main = "Box-Cox and mode relocation")plot(x4, xlab = "sigma", ylab = "xi", ru_scale = TRUE,  main = "Box-Cox, mode relocation and rotation")graphics::par(def_par)

Selecting the Box-Cox parameter in the 1D case

Description

Finds a value of the Box-Cox transformation parameter lambda (\lambda)for which the (positive univariate) random variable with log-density\log f has a density closer to that of a Gaussian randomvariable. Works by estimating a set of quantiles of the distribution impliedby\log f and treating those quantiles as data in a standardBox-Cox analysis. In the following we usetheta (\theta) todenote the argument of\log f on the original scale andphi (\phi) on the Box-Cox transformed scale.

Usage

find_lambda_one_d(  logf,  ...,  ep_bc = 1e-04,  min_phi = ep_bc,  max_phi = 10,  num = 1001,  xdiv = 100,  probs = seq(0.01, 0.99, by = 0.01),  lambda_range = c(-3, 3),  phi_to_theta = NULL,  log_j = NULL)

Arguments

Details

The general idea is to estimate quantiles off correspondingto a set of equally-spaced probabilities inprobs and to use theseestimated quantiles as data in a standard estimation of the Box-Coxtransformation parameter lambda.

The densityf is first evaluated atnum points equallyspaced over the interval (min_phi,max_phi). The continuousdensityf is approximated by attaching trapezium-rule estimates ofprobabilities to the midpoints of the intervals between the points. Afterstandardizing to account for the fact thatf may not be normalized,(min_phi,max_phi) is reset so that values with smallestimated probability (determined byxdiv) are excluded and theprocedure is repeated on this new range. Then the required quantiles areestimated by inferring them from a weighted empirical distributionfunction based on treating the midpoints as data and the estimatedprobabilities at the midpoints as weights.

Value

A list containing the following components

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).

See Also

ru andru_rcpp to performratio-of-uniforms sampling.

find_lambda andfind_lambda_rcppto produce (somewhat) automaticallya list for the argument lambda ofru/ru_rcppfor any value ofd.

find_lambda_one_d_rcpp for a version offind_lambda_one_d that uses the Rcpp package to improveefficiency.

Examples

# Log-normal density ===================# Note: the default value of max_phi = 10 is OK here but this will not# always be the case.lambda <- find_lambda_one_d(logf = dlnorm, log = TRUE)lambdax <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, trans = "BC",        lambda = lambda)# Gamma density ===================alpha <- 1# Choose a sensible value of max_phimax_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.]lambda <- find_lambda_one_d(logf = dgamma, shape = alpha, log = TRUE,                            max_phi = max_phi)lambdax <- ru(logf = dgamma, shape = alpha, log = TRUE, 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(logf = dgamma, shape = alpha, log = TRUE,                            max_phi = max_phi)lambdax <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000,        trans = "BC", lambda = lambda)plot(x)plot(x, ru_scale = TRUE)

Selecting the Box-Cox parameter in the 1D case using Rcpp

Description

Finds a value of the Box-Cox transformation parameter lambda for whichthe (positive univariate) random variable with log-density\log f has a density closer to that of a Gaussian randomvariable. Works by estimating a set of quantiles of the distribution impliedby\log f and treating those quantiles as data in a standardBox-Cox analysis. In the following we usetheta (\theta) todenote the argument of\log f on the original scale andphi (\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

Details

The general idea is to estimate quantiles off correspondingto a set of equally-spaced probabilities inprobs and to use theseestimated quantiles as data in a standard estimation of the Box-Coxtransformation parameterlambda.

The densityf is first evaluated atnum points equally spacedover the interval (min_phi,max_phi). The continuousdensityf is approximated by attaching trapezium-rule estimates ofprobabilities to the midpoints of the intervals between the points. Afterstandardizing to account for the fact thatf may not be normalized,(min_phi,max_phi) is reset so that values with smallestimated probability (determined byxdiv) are excluded and theprocedure is repeated on this new range. Then the required quantiles areestimated by inferring them from a weighted empirical distributionfunction based on treating the midpoints as data and the estimatedprobabilities at the midpoints as weights.

Value

A list containing the following components

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: SeamlessR and C++ Integration.Journal of Statistical Software,40(8), 1-18.doi:10.18637/jss.v040.i08

Eddelbuettel, D. (2013).Seamless R and C++ Integrationwith 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) automaticallya list for the argumentlambda ofru for any value ofd.

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 <- 0sigma <- 1lambda <- find_lambda_one_d_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma)lambdax <- 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_phimax_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)lambdax <- 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)lambdax <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = 1000, trans = "BC",             lambda = lambda)plot(x)plot(x, ru_scale = TRUE)

Selecting the Box-Cox parameter for general d using Rcpp

Description

Finds a value of the Box-Cox transformation parameter lambda for whichthe (positive) random variable with log-density\log f has adensity closer to that of a Gaussian random variable.In the following we usetheta (\theta) to denote the argumentoflogf on the original scale andphi (\phi) on theBox-Cox transformed scale.

Usage

find_lambda_rcpp(  logf,  ...,  d = 1,  n_grid = NULL,  ep_bc = 1e-04,  min_phi = rep(ep_bc, d),  max_phi = rep(10, d),  which_lam = 1:d,  lambda_range = c(-3, 3),  init_lambda = NULL,  phi_to_theta = NULL,  log_j = NULL,  user_args = list())

Arguments

Details

The general idea is to evaluate the densityf on ad-dimensional grid, withn_grid ordinates for each of thed variables.We treat each combination of the variables in the grid as a data pointand perform an estimation of the Box-Cox transformation parameterlambda, in which each data point is weighted by the densityat that point. The vectorsmin_phi andmax_phi define thelimits of the grid andwhich_lam can be used to specify that onlycertain components ofphi are to be transformed.

Value

A list containing the following components

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: SeamlessR and C++ Integration.Journal of Statistical Software,40(8), 1-18.doi:10.18637/jss.v040.i08

Eddelbuettel, D. (2013).Seamless R and C++ Integrationwith Rcpp, Springer, New York. ISBN 978-1-4614-6867-7.

See Also

ru_rcpp to perform ratio-of-uniforms sampling.

find_lambda_one_d_rcpp to produce (somewhat)automatically a list for the argumentlambda ofru for thed = 1 case.

Examples

# Log-normal density ===================# Note: the default value max_phi = 10 is OK here but this will not always# be the caseptr_lnorm <- create_xptr("logdlnorm")mu <- 0sigma <- 1lambda <- find_lambda_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma)lambdax <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, d = 1, n = 1000,             trans = "BC", lambda = lambda)# Gamma density ===================alpha <- 1#  Choose a sensible value of max_phimax_phi <- qgamma(0.999, shape = alpha)# [Of course, typically the quantile function won't be available.  However,# 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_rcpp(logf = ptr_gam, alpha = alpha, max_phi = max_phi)lambdax <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = 1000, trans = "BC",             lambda = lambda)# Generalized Pareto posterior distribution ===================n <- 1000# Sample data from a GP(sigma, xi) distributiongpd_data <- rgpd(m = 100, xi = -0.5, sigma = 1)# Calculate summary statistics for use in the log-likelihoodss <- gpd_sum_stats(gpd_data)# Calculate an initial estimateinit <- c(mean(gpd_data), 0)n <- 1000# Sample on original scale, with no rotation ----------------ptr_gp <- create_xptr("loggp")for_ru_rcpp <- c(list(logf = ptr_gp, init = init, d = 2, n = n,                     lower = c(0, -Inf)), ss, rotate = FALSE)x1 <- do.call(ru_rcpp, for_ru_rcpp)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)# Sample on original scale, with rotation ----------------for_ru_rcpp <- c(list(logf = ptr_gp, init = init, d = 2, n = n,                      lower = c(0, -Inf)), ss)x2 <- do.call(ru_rcpp, for_ru_rcpp)plot(x2, xlab = "sigma", ylab = "xi")abline(a = 0, b = -1 / ss$xm)summary(x2)# Sample on Box-Cox transformed scale ----------------# Find initial estimates for phi = (phi1, phi2),# where phi1 = sigma#   and phi2 = xi + sigma / max(x),# and ranges of phi1 and phi2 over over which to evaluate# the posterior to find a suitable value of lambda.temp <- do.call(gpd_init, ss)min_phi <- pmax(0, temp$init_phi - 2 * temp$se_phi)max_phi <- pmax(0, temp$init_phi + 2 * temp$se_phi)# Set phi_to_theta() that ensures positivity of phi# We use phi1 = sigma and phi2 = xi + sigma / max(data)# Create an external pointer to this C++ functionptr_phi_to_theta_gp <- create_phi_to_theta_xptr("gp")# Note: log_j is set to zero by default inside find_lambda_rcpp()lambda <- find_lambda_rcpp(logf = ptr_gp, ss = ss, d = 2, min_phi = min_phi,                           max_phi = max_phi, user_args = list(xm = ss$xm),                           phi_to_theta = ptr_phi_to_theta_gp)lambda# Sample on Box-Cox transformed, without rotationx3 <- ru_rcpp(logf = ptr_gp, ss = ss, d = 2, n = n, trans = "BC",              lambda = lambda, rotate = FALSE)plot(x3, xlab = "sigma", ylab = "xi")abline(a = 0, b = -1 / ss$xm)summary(x3)# Sample on Box-Cox transformed, with rotationx4 <- ru_rcpp(logf = ptr_gp, ss = ss, d = 2, n = n, trans = "BC",              lambda = lambda)plot(x4, xlab = "sigma", ylab = "xi")abline(a = 0, b = -1 / ss$xm)summary(x4)def_par <- graphics::par(no.readonly = TRUE)par(mfrow = c(2,2), mar = c(4, 4, 1.5, 1))plot(x1, xlab = "sigma", ylab = "xi", ru_scale = TRUE,  main = "mode relocation")plot(x2, xlab = "sigma", ylab = "xi", ru_scale = TRUE,  main = "mode relocation and rotation")plot(x3, xlab = "sigma", ylab = "xi", ru_scale = TRUE,  main = "Box-Cox and mode relocation")plot(x4, xlab = "sigma", ylab = "xi", ru_scale = TRUE,  main = "Box-Cox, mode relocation and rotation")graphics::par(def_par)

Initial estimates for Generalized Pareto parameters

Description

Calculates initial estimates and estimated standard errors (SEs) for thegeneralized Pareto parameters\sigma and\xi based on anassumed random sample from this distribution. Also, calculatesinitial estimates and estimated standard errors for\phi₁ =\sigmaand\phi₂ =\xi + \sigma x_(m), wherex_(m) is the samplemaximum threshold exceedance.

Usage

gpd_init(gpd_data, m, xm, sum_gp = NULL, xi_eq_zero = FALSE, init_ests = NULL)

Arguments

Details

The main aim is to calculate an admissible estimate of\theta, i.e., one at which the log-likelihood is finite (necessaryfor the posterior log-density to be finite) at the estimate, andassociated estimated SEs. These are converted into estimates and SEs for\phi. The latter can be used to set values ofmin_phi andmax_phi for input tofind_lambda.

In the default setting (xi_eq_zero = FALSE andinit_ests = NULL) the methods tried are Maximum LikelihoodEstimation (MLE) (Grimshaw, 1993), Probability-Weighted Moments (PWM)(Hosking and Wallis, 1987) and Linear Combinations of Ratios of Spacings(LRS) (Reiss and Thomas, 2007, page 134) in that order.

For\xi < -1 the likelihood is unbounded, MLE may fail when\xi is not greater than-0.5 and the observed Fisherinformation for(\sigma, \xi) has finite variance only if\xi > -0.25. We use the ML estimate provided thatthe estimate of\xi returned fromgpd_mle is greater than-1. We only use the SE if the MLE of\xi is greater than-0.25.

If either the MLE or the SE are not OK then we try PWM. We use the PWMestimate only if is admissible, and the MLE was not OK. We use the PWM SE,but this will bec(NA, NA) if the PWM estimate of\xi is> 1/2. If the estimate is still not OK then we try LRS. As a lastresort, which will tend to occur only when\xi is strongly negative,we set\xi = -1 and estimate sigma conditional on this.

Value

Ifinit_ests is not supplied by the user, a list is returnedwith components

Ifinit_ests is supplied then only the numeric vectorinit_phi is returned.

References

Grimshaw, S. D. (1993) Computing Maximum Likelihood Estimatesfor the Generalized Pareto Distribution. Technometrics, 35(2), 185-191.and Computing (1991) 1, 129-133.doi:10.1007/BF01889987.

Hosking, J. R. M. and Wallis, J. R. (1987) Parameter and QuantileEstimation for the Generalized Pareto Distribution. Technometrics, 29(3),339-349.doi:10.2307/1269343.

Reiss, R.-D., Thomas, M. (2007) Statistical Analysis of Extreme Valueswith Applications to Insurance, Finance, Hydrology and Other Fields.Birkhauser.doi:10.1007/978-3-7643-7399-3.

See Also

gpd_sum_stats to calculate summary statistics foruse ingpd_loglik.

rgpd for simulation from a generalized Pareto

find_lambda to produce (somewhat) automaticallya list for the argumentlambda ofru.

Examples

# Sample data from a GP(sigma, xi) distributiongpd_data <- rgpd(m = 100, xi = 0, sigma = 1)# Calculate summary statistics for use in the log-likelihoodss <- gpd_sum_stats(gpd_data)# Calculate initial estimatesdo.call(gpd_init, ss)

Generalized Pareto posterior log-density

Description

Calculates the generalized Pareto posterior log-density based on a particularprior for the generalized Pareto parameters, a Maximal Data Information(MDI) prior truncated to\xi \geq -1 in order to produce aposterior density that is proper.

Usage

gpd_logpost(pars, ss)

Arguments

Value

A numeric scalar. The value of the log-likelihood.

References

Northrop, P. J. and Attalides, N. (2016) Posterior propriety inBayesian extreme value analyses using reference priors. Statistica Sinica,26(2), 721-743,doi:10.5705/ss.2014.034.

See Also

gpd_sum_stats to calculate summary statistics foruse ingpd_loglik.

rgpd for simulation from a generalized Pareto

Examples

# Sample data from a GP(sigma, xi) distributiongpd_data <- rgpd(m = 100, xi = 0, sigma = 1)# Calculate summary statistics for use in the log-likelihoodss <- gpd_sum_stats(gpd_data)# Calculate the generalized Pareto log-posteriorgpd_logpost(pars = c(1, 0), ss = ss)

Generalized Pareto summary statistics

Description

Calculates summary statistics involved in the Generalized Paretolog-likelihood.

Usage

gpd_sum_stats(gpd_data)

Arguments

Value

A list with components

See Also

rgpd for simulation from a generalized Paretodistribution.

Examples

# Sample data from a GP(sigma, xi) distributiongpd_data <- rgpd(m = 100, xi = 0, sigma = 1)# Calculate summary statistics for use in the log-likelihoodss <- gpd_sum_stats(gpd_data)

Plot diagnostics for an ru object

Description

plot method for class"ru". Ford = 1 a histogram ofthe simulated values is plotted with a the density function superimposed.The density is normalized crudely using the trapezium rule. Ford = 2 a scatter plot of the simulated values is produced withdensity contours superimposed. Ford > 2 pairwise plots of thesimulated values are produced.

Usage

## S3 method for class 'ru'plot(  x,  y,  ...,  n = ifelse(x$d == 1, 1001, 101),  prob = c(0.1, 0.25, 0.5, 0.75, 0.95, 0.99),  ru_scale = FALSE,  rows = NULL,  xlabs = NULL,  ylabs = NULL,  var_names = NULL,  points_par = list(col = 8))

Arguments

Value

No return value, only the plot is produced.

See Also

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

Examples

# Log-normal density ----------------x <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, lower = 0, init = 1)plot(x)# Improve appearance using arguments to plot() and hist()plot(x, breaks = seq(0, ceiling(max(x$sim_vals)), by = 0.25),  xlim = c(0, 10))# Two-dimensional normal with positive association ----------------rho <- 0.9covmat <- 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)}x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0))plot(x)

Print method for an"ru" object

Description

print method for class"ru".

Usage

## S3 method for class 'ru'print(x, ...)

Arguments

Details

Simply prints the call toru orru_rcpp.

Value

The argumentx, invisibly.

See Also

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

plot.ru for a diagnostic plot.

Generalized Pareto simulation

Description

Simulates a sample of sizem from a generalized Pareto distribution.

Usage

rgpd(m = 1, sigma = 1, xi = 0)

Arguments

Value

A numeric vector. A generalized Pareto sample of sizem.

Examples

# Sample data from a GP(sigma, xi) distributiongpd_data <- rgpd(m = 100, xi = 0, sigma = 1)

Generalized ratio-of-uniforms sampling

Description

Uses the generalized ratio-of-uniforms method to simulate from adistribution with log-density\log f (up to an additiveconstant). The densityf must be bounded, perhaps after atransformation of variable.

Usage

ru(  logf,  ...,  n = 1,  d = 1,  init = NULL,  mode = 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,  shoof = 0.2)

Arguments

Details

For information about the generalised ratio-of-uniforms method andtransformations see theIntroducing rust vignette. This can also be accessed usingvignette("rust-a-vignette", package = "rust").

Iftrans = "none" androtate = FALSE thenruimplements the (multivariate) generalized ratio of uniforms methoddescribed in Wakefield, Gelfand and Smith (1991) using a targetdensity whose mode is relocated to the origin (‘mode relocation’) in thehope of increasing efficiency.

Iftrans = "BC" then marginal Box-Cox transformations of each ofthed variables is performed, with parameters supplied inlambda. The functionphi_to_theta may be used, ifnecessary, to ensure positivity of the variables prior to Box-Coxtransformation.

Iftrans = "user" then the functionphi_to_theta enablesthe user to specify their own transformation.

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

Ifd is greater than one androtate = TRUE then a rotationof the variable axes is performedafter mode relocation. Therotation is based on the Choleski decomposition (seechol) of theestimated Hessian (computed usingoptimHessof the negatedlog-density after any user-supplied transformation or Box-Coxtransformation. If any of the eigenvalues of the estimated Hessian arenon-positive (which may indicate that the estimated mode oflogfis close to a variable boundary) thenrotate is set toFALSEwith a warning. A warning is also given if this happens whend = 1.

The default value of the tuning parameterr is 1/2, which islikely to be close to optimal in many cases, particularly iftrans = "BC".

Value

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

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.doi:10.1007/BF01889987.

See Also

ru_rcpp for a version ofru that usesthe Rcpp package to improve efficiency.

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

plot.ru for a diagnostic plot.

find_lambda_one_d to produce (somewhat) automaticallya list for the argumentlambda ofru for thed = 1 case.

find_lambda to produce (somewhat) automaticallya list for the argumentlambda ofru for any value ofd.

optim for choices of the argumentsa_method,b_method,a_control andb_control.

nlminb for choices of the argumentsa_control andb_control.

optimHess for Hessian estimation.

chol for the Choleski decomposition.

Examples

# 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.9covmat <- 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 <- 0x <- 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 handx <- 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 <- 10x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000,        lower = 0, init = alpha)alpha <- 1x <- 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 <- 1x <- 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 NAphi_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/3x <- 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)# Generalized Pareto posterior distribution ===================# Sample data from a GP(sigma, xi) distributiongpd_data <- rgpd(m = 100, xi = -0.5, sigma = 1)# Calculate summary statistics for use in the log-likelihoodss <- gpd_sum_stats(gpd_data)# Calculate an initial estimateinit <- c(mean(gpd_data), 0)# Mode relocation only ----------------n <- 1000x1 <- 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$paplot(x, ru_scale = TRUE)# Example 4 from Wakefield et al. (1991) ===================# Bivariate normal x bivariate student-tlog_norm_t <- function(x, mean = rep(0, d), sigma1 = diag(d), sigma2 = diag(d)) {  x <- matrix(x, ncol = length(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.9covmat <- matrix(c(1, rho, rho, 1), 2, 2)y <- c(0, 0)# Case in the top right corner of Table 3x <- 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 acceptancex <- ru(logf = log_norm_t, mean = y, sigma1 = covmat, sigma2 = covmat,  d = 2, n = 10000, init = y, rotate = TRUE)x$pa# Normal x log-normal: different Box-Cox parameters ==================norm_lognorm <- function(x, ...) {  dnorm(x[1], ...) + dlnorm(x[2], ...)}x <- ru(logf = norm_lognorm, log = TRUE, n = 1000, d = 2, init = c(-1, 0),        trans = "BC", lambda = c(1, 0))plot(x)plot(x, ru_scale = TRUE)

Generalized ratio-of-uniforms sampling using C++ via Rcpp

Description

Uses the generalized ratio-of-uniforms method to simulate from adistribution with log-density\log f (up to an additiveconstant). The densityf must be bounded, perhaps after atransformation of variable.The fileuser_fns.cpp that is sourced before running the examplesbelow is available at the rust Github page athttps://raw.githubusercontent.com/paulnorthrop/rust/master/src/user_fns.cpp.

Usage

ru_rcpp(  logf,  ...,  n = 1,  d = 1,  init = NULL,  mode = 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,  shoof = 0.2)

Arguments

Details

For information about the generalised ratio-of-uniforms method andtransformations see theIntroducing rust vignette. See alsoRusting faster: Simulation using Rcpp.

These vignettes can also be accessed usingvignette("rust-a-vignette", package = "rust") andvignette("rust-c-using-rcpp-vignette", package = "rust").

Iftrans = "none" androtate = FALSE thenruimplements the (multivariate) generalized ratio of uniforms methoddescribed in Wakefield, Gelfand and Smith (1991) using a targetdensity whose mode is relocated to the origin (‘mode relocation’) in thehope of increasing efficiency.

Iftrans = "BC" then marginal Box-Cox transformations of each ofthed variables is performed, with parameters supplied inlambda. The functionphi_to_theta may be used, ifnecessary, to ensure positivity of the variables prior to Box-Coxtransformation.

Iftrans = "user" then the functionphi_to_theta enablesthe user to specify their own transformation.

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

Ifd is greater than one androtate = TRUE then a rotationof the variable axes is performedafter mode relocation. Therotation is based on the Choleski decomposition (seechol) of theestimated Hessian (computed usingoptimHessof the negatedlog-density after any user-supplied transformation or Box-Coxtransformation. If any of the eigenvalues of the estimated Hessian arenon-positive (which may indicate that the estimated mode oflogfis close to a variable boundary) thenrotate is set toFALSEwith a warning. A warning is also given if this happens whend = 1.

The default value of the tuning parameterr is 1/2, which islikely to be close to optimal in many cases, particularly iftrans = "BC".

Value

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

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.doi:10.1007/BF01889987.

Eddelbuettel, D. and Francois, R. (2011). Rcpp: SeamlessR and C++ Integration.Journal of Statistical Software,40(8), 1-18.doi:10.18637/jss.v040.i08

Eddelbuettel, D. (2013).Seamless R and C++ Integrationwith Rcpp, Springer, New York. ISBN 978-1-4614-6867-7.

See Also

ru for a version ofru_rcpp thataccepts R functions as arguments.

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

plot.ru for a diagnostic plot.

find_lambda_one_d_rcpp to produce (somewhat)automatically a list for the argumentlambda ofru for thed = 1 case.

find_lambda_rcpp to produce (somewhat) automaticallya list for the argumentlambda ofru for any value ofd.

optim for choices of the argumentsa_method,b_method,a_control andb_control.

nlminb for choices of the argumentsa_control andb_control.

optimHess for Hessian estimation.

chol for the Choleski decomposition.

Examples

n <- 1000# Normal density ===================# One-dimensional standard normal ----------------ptr_N01 <- create_xptr("logdN01")x <- ru_rcpp(logf = ptr_N01, d = 1, n = n, init = 0.1)# Two-dimensional standard normal ----------------ptr_bvn <- create_xptr("logdnorm2")rho <- 0x <- ru_rcpp(logf = ptr_bvn, rho = rho, d = 2, n = n,  init = c(0, 0))# Two-dimensional normal with positive association ===================rho <- 0.9# No rotation.x <- ru_rcpp(logf = ptr_bvn, rho = rho, d = 2, n = n, init = c(0, 0),             rotate = FALSE)# With rotation.x <- ru_rcpp(logf = ptr_bvn, rho = rho, d = 2, n = n, init = c(0, 0))# Using general multivariate normal function.ptr_mvn <- create_xptr("logdmvnorm")covmat <- matrix(rho, 2, 2) + diag(1 - rho, 2)x <- ru_rcpp(logf = ptr_mvn, sigma = covmat, d = 2, n = n, init = c(0, 0))# Three-dimensional normal with positive association ----------------covmat <- matrix(rho, 3, 3) + diag(1 - rho, 3)# No rotation.x <- ru_rcpp(logf = ptr_mvn, sigma = covmat, d = 3, n = n,             init = c(0, 0, 0), rotate = FALSE)# With rotation.x <- ru_rcpp(logf = ptr_mvn, sigma = covmat, d = 3, n = n,             init = c(0, 0, 0))# Log-normal density ===================ptr_lnorm <- create_xptr("logdlnorm")mu <- 0sigma <- 1# Sampling on original scale ----------------x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, d = 1, n = n,             lower = 0, init = exp(mu))# Box-Cox transform with lambda = 0 ----------------lambda <- 0x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, d = 1, n = n,             lower = 0, init = exp(mu), trans = "BC", lambda = lambda)# Equivalently, we could use trans = "user" and supply the (inverse) Box-Cox# transformation and the log-Jacobian by handptr_phi_to_theta_lnorm <- create_phi_to_theta_xptr("exponential")ptr_log_j_lnorm <- create_log_j_xptr("neglog")x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, d = 1, n = n,  init = 0.1, trans = "user", phi_to_theta = ptr_phi_to_theta_lnorm,  log_j = ptr_log_j_lnorm)# 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 ----------------ptr_gam <- create_xptr("logdgamma")alpha <- 10x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = n,  lower = 0, init = alpha)alpha <- 1x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = n,  lower = 0, init = alpha)# Box-Cox transform with lambda = 1/3 works well for shape >= 1. -----------alpha <- 1x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = n,  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 handlambda <- 1/3ptr_phi_to_theta_bc <- create_phi_to_theta_xptr("bc")ptr_log_j_bc <- create_log_j_xptr("bc")x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = n,  trans = "user", phi_to_theta = ptr_phi_to_theta_bc, log_j = ptr_log_j_bc,  user_args = list(lambda = lambda), init = alpha)summary(x)# Generalized Pareto posterior distribution ===================# Sample data from a GP(sigma, xi) distributiongpd_data <- rgpd(m = 100, xi = -0.5, sigma = 1)# Calculate summary statistics for use in the log-likelihoodss <- gpd_sum_stats(gpd_data)# Calculate an initial estimateinit <- c(mean(gpd_data), 0)n <- 1000# Mode relocation only ----------------ptr_gp <- create_xptr("loggp")for_ru_rcpp <- c(list(logf = ptr_gp, init = init, d = 2, n = n,                 lower = c(0, -Inf)), ss, rotate = FALSE)x1 <- do.call(ru_rcpp, for_ru_rcpp)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 ----------------for_ru_rcpp <- c(list(logf = ptr_gp, init = init, d = 2, n = n,                 lower = c(0, -Inf)), ss)x2 <- do.call(ru_rcpp, for_ru_rcpp)plot(x2, xlab = "sigma", ylab = "xi")abline(a = 0, b = -1 / ss$xm)summary(x2)# Cauchy ========================ptr_c <- create_xptr("logcauchy")# 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_rcpp(logf = ptr_c, log = TRUE, init = 0, r = 1.26, n = 1000)# Half-Cauchy ===================ptr_hc <- create_xptr("loghalfcauchy")# 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_rcpp(logf = ptr_hc, init = 0, trans = "BC", lambda = 0, n = 1000)x$paplot(x, ru_scale = TRUE)# Example 4 from Wakefield et al. (1991) ===================# Bivariate normal x bivariate student-tptr_normt <- create_xptr("lognormt")rho <- 0.9covmat <- matrix(c(1, rho, rho, 1), 2, 2)y <- c(0, 0)# Case in the top right corner of Table 3x <- ru_rcpp(logf = ptr_normt, mean = y, sigma1 = covmat, sigma2 = covmat,  d = 2, n = 10000, init = y, rotate = FALSE)x$pa# Rotation increases the probability of acceptancex <- ru_rcpp(logf = ptr_normt, mean = y, sigma1 = covmat, sigma2 = covmat,  d = 2, n = 10000, init = y, rotate = TRUE)x$pa

Internal rust functions

Description

Internal rust functions

Usage

box_cox(x, lambda = 1, gm = 1, lambda_tol = 1e-06)box_cox_vec(x, lambda = 1, gm = 1, lambda_tol = 1e-06)optim_box_cox(  x,  w,  lambda_range = c(-3, 3),  start = NULL,  which_lam = 1:ncol(x))n_grid_fn(d)init_psi_calc(phi_to_psi, phi, lambda, gm, w, which_lam)log_gpd_mdi_prior(pars)gpd_loglik(pars, gpd_data, m, xm, sum_gp)gpd_mle(gpd_data)fallback_gp_mle(init, ...)gpd_obs_info(gpd_pars, y, eps = 1e-05, m = 3)find_a(  neg_logf_rho,  init_psi,  d,  r,  lower,  upper,  algor,  method,  control,  shoof,  ...)find_bs(  f_rho,  d,  r,  lower,  upper,  f_mode,  ep,  vals,  conv,  algor,  method,  control,  shoof,  ...)cpp_find_a(  init_psi,  lower,  upper,  algor,  method,  control,  a_obj_fun,  ru_args,  shoof)cpp_find_bs(  lower,  upper,  ep,  vals,  conv,  algor,  method,  control,  lower_box_fun,  upper_box_fun,  ru_args,  shoof)fac3(d)

Details

These functions are not intended to be called by the user.

Summarizing ratio-of-uniforms samples

Description

summary method for class"ru".

print method for an objectobject of class"summary.ru".

Usage

## S3 method for class 'ru'summary(object, ...)## S3 method for class 'summary.ru'print(x, ...)

Arguments

Value

Forsummary.lm: a list of the following components fromobject:

• information about the ratio-of-uniforms bounding box, i.e.,object$box

• an estimate of the probability of acceptance, i.e.,object$pa

• a summary of the simulated values, viasummary(object$sim_vals)

Forprint.summary.ru: the argumentx, invisibly.

See Also

ru for descriptions ofobject$sim_vals andobject$box.

plot.ru for a diagnostic plot.

Examples

# one-dimensional standard normal ----------------x <- ru(logf = function(x) -x ^ 2 / 2, d = 1, n = 1000, init = 0)summary(x)# two-dimensional normal with positive association ----------------rho <- 0.9covmat <- 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)}x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0))summary(x)