# The right hand side of the log likelihood of a Batschelet-type distribution
ll_rhs_bat <- function(x, mu, kp, lam, tlam_fun) {
kp * sum(cos(tlam_fun(x - mu, lam)))
}
# # The left hand side of the log likelihood of a inverse Batschelet distribution
# ll_lhs_invbat <- function(n, kp, lam) {
# -n * (logBesselI(kp, 0) + log(K_kplam(kp, lam)))
# }
sample_mu_bat <- function(x, mu_cur, kp, lam, tlam_fun, mu_logprior_fun) {
# Sample a candidate from the distribution of mu when we have a von Mises
# distribution.
C_j <- sum(cos(x))
S_j <- sum(sin(x))
R_j <- sqrt(C_j^2 + S_j^2)
mu_can <- circglmbayes::rvmc(1, mu_cur, R_j * kp)
ll_can <- ll_rhs_bat(x, mu_can, kp, lam, tlam_fun)
ll_cur <- ll_rhs_bat(x, mu_cur, kp, lam, tlam_fun)
# The proposal is von Mises and thus symmetric, so the transition
# probabilities of MH are omitted here.
mu_lograt <- ll_can + mu_logprior_fun(mu_can) - ll_cur - mu_logprior_fun(mu_cur)
if (mu_lograt > log(stats::runif(1))) {
return(mu_can)
} else {
return(mu_cur)
}
}
sample_mu_bat_2 <- function(x, mu_cur, kp, lam, tlam_fun, mu_logprior_fun) {
# Sample a candidate from the distribution of mu when we have a von Mises
# distribution.
C_j <- sum(cos(x))
S_j <- sum(sin(x))
R_j <- sqrt(C_j^2 + S_j^2)
mu_hat <- atan2(S_j, C_j)
mu_can <- circglmbayes::rvmc(1, mu_hat, R_j * kp)
ll_can <- ll_rhs_bat(x, mu_can, kp, lam, tlam_fun)
ll_cur <- ll_rhs_bat(x, mu_cur, kp, lam, tlam_fun)
logp_mu_can_to_cur <- dvm(mu_cur, mu_hat, kp, log = TRUE)
logp_mu_cur_to_can <- dvm(mu_can, mu_hat, kp, log = TRUE)
mu_lograt <- ll_can + mu_logprior_fun(mu_can) + logp_mu_can_to_cur -
ll_cur - mu_logprior_fun(mu_cur) - logp_mu_cur_to_can
if (mu_lograt > log(stats::runif(1))) {
return(mu_can)
} else {
return(mu_cur)
}
}
# Reparametrized gamma proposal to make tuning parameter and mean interpretable.
# This is equal to chi square with 'mean' degrees of freedom if var_tune = 1.
dgammaprop <- function(x, mean = 1, var_tune = 1, log = FALSE) {
gamma_var <- 2 * mean * var_tune
gamma_scale <- gamma_var / mean
gamma_shape <- mean / gamma_scale
stats::dgamma(x, shape = gamma_shape, scale = gamma_scale, log = log)
}
rgammaprop <- function(n = 1, mean = 1, var_tune = 1) {
gamma_var <- 2 * mean * var_tune
gamma_scale <- gamma_var / mean
gamma_shape <- mean / gamma_scale
stats::rgamma(n = n, shape = gamma_shape, scale = gamma_scale)
}
sample_kp_bat <- function(x, mu, kp_cur, lam, llbat, kp_logprior_fun, var_tune = 1) {
# Sample a candidate
if (kp_cur > 0) {
kp_can <- rgammaprop(1, mean = kp_cur, var_tune = var_tune)
logp_kp_can_to_cur <- dgammaprop(kp_cur, kp_can, var_tune, log = TRUE)
logp_kp_cur_to_can <- dgammaprop(kp_can, kp_cur, var_tune, log = TRUE)
} else {
# If kp_cur == 0, usual sampling from gamma breaks down, so we retry by
# drawing proposals from the distribution with kp_cur == .1.
kp_can <- rgammaprop(1, mean = .1, var_tune = var_tune)
logp_kp_can_to_cur <- dgammaprop(.1, kp_can, var_tune, log = TRUE)
logp_kp_cur_to_can <- dgammaprop(kp_can, .1, var_tune, log = TRUE)
}
ll_can <- llbat(x, mu, kp_can, lam, log = TRUE)
ll_cur <- llbat(x, mu, kp_cur, lam, log = TRUE)
kp_lograt <- ll_can + kp_logprior_fun(kp_can) + logp_kp_can_to_cur -
ll_cur - kp_logprior_fun(kp_cur) - logp_kp_cur_to_can
if (kp_lograt > log(stats::runif(1))) {
return(kp_can)
} else {
return(kp_cur)
}
}
sample_lam_bat <- function(x, mu, kp, lam_cur, llbat, lam_logprior_fun, lam_bw = .05) {
# Sample a candidate
lam_can <- stats::runif(1, max(-1, lam_cur - lam_bw), min(1, lam_cur + lam_bw))
ll_can <- llbat(x, mu, kp, lam_can)
ll_cur <- llbat(x, mu, kp, lam_cur)
logp_lam_can_to_cur <- stats::dunif(lam_cur, max(-1, lam_can - lam_bw),
min(1, lam_can + lam_bw), log = TRUE)
logp_lam_cur_to_can <- stats::dunif(lam_can, max(-1, lam_cur - lam_bw),
min(1, lam_cur + lam_bw), log = TRUE)
lam_lograt <- ll_can + lam_logprior_fun(lam_can) + logp_lam_can_to_cur -
ll_cur - lam_logprior_fun(lam_cur) - logp_lam_cur_to_can
if (lam_lograt > log(stats::runif(1))) {
return(lam_can)
} else {
return(lam_cur)
}
}
sample_kp_and_lam_bat <- function(x, mu, kp_cur, lam_cur, llbat, lam_bw = .05,
kp_logprior_fun, lam_logprior_fun, var_tune = 1) {
if (kp_cur > 0) {
kp_can <- rgammaprop(1, mean = kp_cur, var_tune = var_tune)
logp_kp_can_to_cur <- dgammaprop(kp_cur, kp_can, var_tune, log = TRUE)
logp_kp_cur_to_can <- dgammaprop(kp_can, kp_cur, var_tune, log = TRUE)
} else {
# If kp_cur == 0, usual sampling from gamma breaks down, so we retry by
# drawing proposals from the distribution with kp_cur == .1.
kp_can <- rgammaprop(1, mean = .1, var_tune = var_tune)
logp_kp_can_to_cur <- dgammaprop(.1, kp_can, var_tune, log = TRUE)
logp_kp_cur_to_can <- dgammaprop(kp_can, .1, var_tune, log = TRUE)
}
lam_can <- stats::runif(1,
max(-1, lam_cur - lam_bw),
min(1, lam_cur + lam_bw))
logp_kp_can_to_cur <- dgammaprop(kp_cur, kp_can, var_tune, log = TRUE)
logp_kp_cur_to_can <- dgammaprop(kp_can, kp_cur, var_tune, log = TRUE)
logp_lam_can_to_cur <- stats::dunif(lam_cur,
max(-1, lam_can - lam_bw),
min(1, lam_can + lam_bw),
log = TRUE)
logp_lam_cur_to_can <- stats::dunif(lam_can,
max(-1, lam_cur - lam_bw),
min(1, lam_cur + lam_bw),
log = TRUE)
ll_can <- llbat(x, mu, kp_can, lam_can, log = TRUE)
ll_cur <- llbat(x, mu, kp_cur, lam_cur, log = TRUE)
kplam_lograt <- ll_can + kp_logprior_fun(kp_can) + lam_logprior_fun(lam_can) +
logp_kp_can_to_cur + logp_lam_can_to_cur -
ll_cur - kp_logprior_fun(kp_cur) - lam_logprior_fun(lam_cur) -
logp_kp_cur_to_can - logp_lam_cur_to_can
if (kplam_lograt > log(stats::runif(1))) {
return(c(kp_can, lam_can))
} else {
return(c(kp_cur, lam_cur))
}
}
#' A rescaled beta prior for lambda
#'
#' This prior is symmetric between -1 and 1, and captures the prior belief that
#' values of \code{lam} on the boundary of the parameter space are a prior
#' unlikely.
#'
#' @param lam Numeric;
#'
#' @return The log-prior probability.
#' @export
#'
lam_beta_log_prior_2_2 <- function(lam) {
stats::dbeta( (lam + 1) / 2, 2, 2, log = TRUE)
}
logBesselI <- function(x, nu) log(besselI(x, nu, expon.scaled = TRUE)) + x
#' The Jeffreys prior for the von Mises distribution
#'
#' @param kp Numeric;
#'
#' @return An unnormalized value.
#' @export
vm_kp_jeffreys_prior <- function(kp) {
logAkp <- logBesselI(kp, 1) - logBesselI(kp, 0)
return(
exp(0.5 * (log(kp) +
logAkp +
log(0.5 +
exp(logBesselI(kp, 2) - log(2) - logBesselI(kp, 0)) -
exp(logAkp) ^ 2))))
}
#' @rdname vm_kp_jeffreys_prior
#' @export
vm_kp_jeffreys_logprior <- function(kp) {
logAkp <- logBesselI(kp, 1) - logBesselI(kp, 0)
return(
0.5 * (log(kp) +
logAkp +
log(0.5 +
exp(logBesselI(kp, 2) - log(2) - logBesselI(kp, 0)) -
exp(logAkp) ^ 2)))
}
#' MCMC sampling for Batschelet-type distributions.
#'
#' @param x A numeric vector of angles, in radians
#' @param Q Integer; The number of iterations to return after taking burn in and
#' thinning into account.
#' @param burnin Integer; The number of (non-thinned) iterations to discard. No
#' burn in is performed by default.
#' @param thin Integer; Number of iterations to sample for each saved iteration.
#' Defaults to 1, which means no thinning.
#' @param n_comp Integer; Fixed number of components to estimate.
#' @param bat_type Either 'inverse' or 'power', the type of distribution to fit.
#' The two distributions are similar, but the power Batschelet distribution is
#' computationally much less demanding.
#' @param init_pmat A numeric matrix with \code{n_comp} rows and four columns,
#' corresponding to \eqn{\mu, \kappa, \lambda, \alpha}, in that order. Gives
#' starting values for the parameters. If any element is \code{NA}, it will be
#' given a default starting value. For \eqn{mu}, the default starting values
#' are equally spaced on the circle. For \eqn{\kappa}, the default starting
#' value is 5. For \eqn{\lambda}, the default starting value is 0, which
#' corresponds to the von Mises distribution. For \eqn{\alpha}, the default
#' starting value is \code{1/n_comp}.
#' @param fixed_pmat A numeric matrix with \code{n_comp} rows and four columns,
#' corresponding to \eqn{\mu, \kappa, \lambda, \alpha}, in that order. Any
#' element that is not \code{NA} in this matrix will be held constant at the
#' given value and not sampled.
#' @param mu_logprior_fun Function; A function with a single argument, which
#' returns the log of the prior probability of \eqn{\mu}. Defaults to a
#' uniform prior function.
#' @param kp_logprior_fun Function; A function with a single argument, which
#' returns the log of the prior probability of \eqn{\kappa}. Defaults to a
#' uniform prior function. In contrast to the other parameters, for
#' \eqn{\kappa} the constant (uniform) prior is improper.
#' @param lam_logprior_fun Function; A function with a single argument, which
#' returns the log of the prior probability of \eqn{\lambda}. Defaults to a
#' uniform prior function.
#' @param alph_prior_param Integer vector; The mixture weight parameter vector
#' \eqn{\alpha} is given its conjugate Dirichlet prior. The default is
#' \code{rep(1, n_comp)}, which is the noninformative uniform prior over the
#' \code{n_comp} simplex.
#' @param joint_kp_lam Logical; If \code{TRUE}, the parameters \code{kp} and
#' \code{lam} are drawn jointly. This can be beneficial if these are strongly
#' correlated.
#' @param verbose Integer up to 4; Determines the amount of printed debug
#' information.
#' @param lam_bw Numeric; the maximum distance from the current lambda at which
#' uniform proposals are drawn.
#' @param kp_bw Numeric; A tuning parameter for kappa proposals. If \code{kp_bw
#' == 1}, the chi-square distribution is used. Often, this distribution is too
#' wide, so this parameter can be set to \code{0 < kp_bw < 1} to use a gamma
#' proposal which has lower variance than the chi-square.
#' @param compute_variance Logical; Whether to add circular variance to the
#' returned mcmc sample.
#' @param compute_waic Logical; Whether to compute the WAIC. Can be
#' computationally demanding if \code{n * Q} is large.
#'
#' @importFrom stats var
#'
#' @return A numeric matrix of sampled parameter values.
#' @export
#'
#' @examples
#' x <- rinvbatmix(100)
#' mcmcBatscheletMixture(x, Q = 10)
mcmcBatscheletMixture <- function(x, Q = 1000,
burnin = 0, thin = 1,
n_comp = 4,
bat_type = "inverse",
init_pmat = matrix(NA, n_comp, 4),
fixed_pmat = matrix(NA, n_comp, 4),
joint_kp_lam = FALSE,
kp_bw = 1,
lam_bw = .05,
mu_logprior_fun = function(mu) -log(2*pi),
kp_logprior_fun = function(kp) 1,
lam_logprior_fun = function(lam) -log(2),
alph_prior_param = rep(1, n_comp),
compute_variance = TRUE,
compute_waic = FALSE,
verbose = 0) {
# Select Batschelet type
if (bat_type == "inverse") {
dbat_fun <- dinvbat
llbat <- likinvbat
tlam_fun <- t_lam
} else if (bat_type == "power") {
dbat_fun <- dpowbat
llbat <- likpowbat
tlam_fun <- tpow_lam
} else {
stop("Unknown Batschelet type.")
}
# A matrix with logicals for each initial value of the parameter matrix.
na_fixedpmat <- is.na(fixed_pmat)
na_initpmat <- is.na(init_pmat)
if (any(!na_fixedpmat[, 2] & fixed_pmat[,2] < 0)) stop("Invalid fixed kappa value.")
if (any(!na_fixedpmat[, 3] & abs(fixed_pmat[,3]) > 1)) stop("Invalid fixed lambda value.")
# Set initial values if the initial parameter matrix is not given for that parameter (has NAs).
init_pmat[, 1] <- ifelse(na_initpmat[, 1], seq(0, 2*pi, length.out = n_comp + 1)[-1], init_pmat[, 1])
init_pmat[, 2] <- ifelse(na_initpmat[, 2], rep(5, n_comp), init_pmat[, 2])
init_pmat[, 3] <- ifelse(na_initpmat[, 3], rep(0, n_comp), init_pmat[, 3])
init_pmat[, 4] <- ifelse(na_initpmat[, 4], rep(1/n_comp, n_comp), init_pmat[, 4])
# Force x to be in range -pi, pi.
x <- force_neg_pi_pi(x)
n <- length(x)
# Initialize parameters
mu_cur <- init_pmat[, 1]
kp_cur <- init_pmat[, 2]
lam_cur <- init_pmat[, 3]
alph_cur <- init_pmat[, 4]
# Matrix of acceptance totals for kp and lam, which will be divided by the
# total later.
acc_mat <- matrix(0, nrow = n_comp, ncol = 2)
colnames(acc_mat) <- c("kp", "lam")
rownames(acc_mat) <- 1:n_comp
# Initialize latent group labeling
z_cur <- integer(n)
output_matrix <- matrix(NA, nrow = Q, ncol = n_comp*4)
colnames(output_matrix) <- c(paste0("mu_", 1:n_comp),
paste0("kp_", 1:n_comp),
paste0("lam_", 1:n_comp),
paste0("alph_", 1:n_comp))
ll_vec <- numeric(Q)
# Add WAIC log-likelihood accumulation vector.
if (compute_waic) {
ll_each_th_curpars <- matrix(NA, nrow = Q, ncol = n)
}
if (compute_variance) {
variance_matrix <- matrix(NA, nrow = Q, ncol = n_comp*3)
colnames(variance_matrix) <- c(paste0("mean_res_len_", 1:n_comp),
paste0("circ_var_", 1:n_comp),
paste0("circ_sd_", 1:n_comp))
}
Qbythin <- Q * thin + burnin
if (verbose) cat("Starting MCMC sampling.\n")
if (verbose > 1) cat("Iteration:\n")
for (i in 1:Qbythin) {
if (verbose > 1) cat(sprintf("%5s, ", i))
### Sample group assignments z
# Probability of each component for each data point.
W <- sapply(1:n_comp, function(k) {
alph_cur[k] * dbat_fun(x, mu_cur[k], kp_cur[k], lam_cur[k])
})
w_rowsum <- rowSums(W)
# If some datapoints are numerically equal to zero, assign it equally to
# each component. Hopefully, this does not happen again in the next
# iteration.
W[w_rowsum == 0, ] <- 1/n_comp
w_rowsum[w_rowsum == 0] <- 1
W <- W / w_rowsum
# Randomly sample group assignments from the component probabilities.
z_cur <- apply(W, 1, function(row_probs) sample(x = 1:n_comp, size = 1,
prob = row_probs))
# Sample weights alph
dir_asum <- sapply(1:n_comp, function(j) sum(z_cur == j))
alph_cur <- ifelse(na_fixedpmat[, 4],
MCMCpack::rdirichlet(1, dir_asum + alph_prior_param),
alph_cur)
# After sampling the current group assignments, the parameters for each
# component only can be sampled separately.
for (j in 1:n_comp) {
if (verbose > 2) cat(j)
# Dataset assigned to this component.
x_j <- x[z_cur == j]
# Check whether anything is assigned to this components. If not, don't
# update the parameters.
if (length(x_j) == 0) {
if (verbose > 1) cat("---")
next
}
if (verbose > 2) cat("m")
# Sample mu
if (na_fixedpmat[j, 1]) {
mu_cur[j] <- sample_mu_bat_2(x_j, mu_cur[j], kp_cur[j], lam_cur[j],
tlam_fun, mu_logprior_fun)
}
if (joint_kp_lam) {
if (verbose > 2) cat("kl")
kplam_curj <- sample_kp_and_lam_bat(x_j,
mu_cur[j], kp_cur[j], lam_cur[j],
llbat, lam_bw = lam_bw,
kp_logprior_fun, lam_logprior_fun,
var_tune = kp_bw)
# Assign the new values if they are new, and set acceptance count.
if (kp_cur[j] != kplam_curj[1]) {
kp_cur[j] <- kplam_curj[1]
acc_mat[j, 1] <- acc_mat[j, 1] + 1
}
if (lam_cur[j] != kplam_curj[2]) {
lam_cur[j] <- kplam_curj[2]
acc_mat[j, 2] <- acc_mat[j, 2] + 1
}
} else {
if (verbose > 2) cat("k")
# Sample kp
if (na_fixedpmat[j, 2]) {
kp_new <- sample_kp_bat(x_j, mu_cur[j], kp_cur[j], lam_cur[j],
llbat, kp_logprior_fun, var_tune = kp_bw)
if (kp_cur[j] != kp_new) {
kp_cur[j] <- kp_new
acc_mat[j, 1] <- acc_mat[j, 1] + 1
}
}
if (verbose > 2) cat("l")
# Sample lam
if (na_fixedpmat[j, 3]) {
lam_new <- sample_lam_bat(x_j, mu_cur[j], kp_cur[j], lam_cur[j],
llbat, lam_logprior_fun, lam_bw = lam_bw)
if (lam_cur[j] != lam_new) {
lam_cur[j] <- lam_new
acc_mat[j, 2] <- acc_mat[j, 2] + 1
}
}
}
}
# Possibly extremely detailed debugging information.
if (verbose > 3) {
cat("\n",
sprintf("mu: %8s, ", round(mu_cur, 3)), "\n",
sprintf("kp: %8s, ", round(kp_cur, 3)), "\n",
sprintf("lam: %8s, ", round(lam_cur, 3)), "\n",
sprintf("alph: %8s, ", round(alph_cur, 3)), "\n")
}
if (i %% 50 == 0 && verbose == 1) cat(i, ", \n", sep = "")
if (i %% 5 == 0 && verbose > 1) cat("\n")
if (i %% thin == 0 && i >= burnin) {
isav <- (i - burnin) / thin
output_matrix[isav, ] <- c(mu_cur, kp_cur, lam_cur, alph_cur)
# A vector with log-densities for each data point.
each_th_ll <- dbatmix(x = x, dbat_fun = dbat_fun,
mu_cur, kp_cur, lam_cur, alph_cur,
log = TRUE)
ll_vec[isav] <- sum(each_th_ll)
if (compute_waic) {
ll_each_th_curpars[isav, ] <- each_th_ll
}
if (compute_variance) {
# Compute mean resultant lengths for each component.
R_bar_cur <- sapply(1:n_comp, function(i) {
computeMeanResultantLengthBat(kp_cur[i], lam_cur[i], bat_type)
})
# Compute variances.
circ_var_cur <- 1 - R_bar_cur
circ_sd_cur <- computeCircSD(R_bar_cur)
# Put the results in an output matrix.
variance_matrix[isav, ] <- c(R_bar_cur, circ_var_cur, circ_sd_cur)
}
}
}
# Compute acceptance ratio.
acc_mat <- acc_mat / Q
if (verbose) cat("\nFinished.\n")
# Collect output
if (compute_variance) output_matrix <- cbind(output_matrix, variance_matrix)
# The log-posterior function that we have just sampled from.
log_posterior <- function(pvec, data = x) {
# If there is one value in pvec missing, assume it is one of the alpha
# weights because this might happen when bridgesampling for example.
if ((length(pvec) %% 4) == 3) {
n_comp <- (length(pvec) + 1)/4
pvec <- c(pvec, 1 - sum(pvec[(3*n_comp + 1):(4*n_comp - 1)]))
}
n_comp <- length(pvec)/4
mus <- pvec[1:n_comp]
kps <- pvec[(n_comp + 1):(2*n_comp)]
lams <- pvec[(2*n_comp + 1):(3*n_comp)]
alphs <- pvec[(3*n_comp + 1):(4*n_comp)]
if (!identical(sum(alphs), 1)) {
warning("Log-posterior adapting weight vector alpha to sum to one.")
alphs <- alphs / sum(alphs)
}
ll_part <- sum(dbatmix(x, dbat_fun = dbat_fun,
mus, kps, lams, alphs,
log = TRUE))
prior_part <- sum(c(vapply(mus, mu_logprior_fun, 0),
vapply(kps, kp_logprior_fun, 0),
vapply(lams, lam_logprior_fun, 0),
log(MCMCpack::ddirichlet(alphs,
alpha = alph_prior_param))))
ll_part + prior_part
}
# Create a new environment for log_posterior so that the file size of the
# resulting object is not too large.
log_post_env <- new.env()
log_post_env$data <- x
log_post_env$n_comp <- n_comp
log_post_env$dbat_fun <- dbat_fun
log_post_env$mu_logprior_fun <- mu_logprior_fun
log_post_env$kp_logprior_fun <- kp_logprior_fun
log_post_env$lam_logprior_fun <- lam_logprior_fun
log_post_env$alph_prior_param <- alph_prior_param
environment(log_posterior) <- log_post_env
out_list <- list(mcmc_sample = coda::mcmc(output_matrix,
start = burnin + 1,
end = Qbythin,
thin = thin),
ll_vec = ll_vec,
log_posterior = log_posterior,
acceptance_rates = acc_mat,
prior_list = list(mu_logprior_fun, kp_logprior_fun,
lam_logprior_fun, alph_prior_param))
if (compute_waic) {
lppd <- sum(log(colSums(exp(ll_each_th_curpars)))) - n * log(Q)
waic_logofmean <- log(colMeans(exp(ll_each_th_curpars)))
waic_meanoflog <- colMeans(ll_each_th_curpars)
p_waic1 <- 2 * sum(waic_logofmean - waic_meanoflog)
p_waic2 <- sum(apply(ll_each_th_curpars, 2, var))
waic1 <- -2 * (lppd - p_waic1)
waic2 <- -2 * (lppd - p_waic2)
out_list$ic <- list(lppd = lppd,
waic_1 = c(p_waic1 = 2 * p_waic1, waic1 = waic1),
waic_2 = c(p_waic2 = 2 * p_waic2, waic2 = waic2))
} else {
out_list$ic <- list()
}
return(out_list)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.