# Log-normal density ===================
# Note: the default value max_phi = 10 is OK here but this will not always
# be the case
lambda <- find_lambda(logf = dlnorm, log = TRUE)
lambda
x <- ru(logf = dlnorm, 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)
# [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)
lambda
x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000,
trans = "BC", lambda = lambda)
# \donttest{
# 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)
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) 0
lambda <- 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 rotation
x3 <- 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 rotation
x4 <- 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)
# }
Run the code above in your browser using DataLab