Nothing
# find_lscale <- function(x) {
# y <- x
# y[1:2] <- log(x[1:2])
# y
# }
find_expscale <- function(x) {
y <- x
y[1:2] <- exp(x[1:2])
y
}
find_expscale_uni <- function(x) {
y <- x
y[1] <- exp(x[1])
y
}
find_lscale_mat <- function(x) {
y <- x
y[1:2, ] <- log(x[1:2, ])
y
}
find_lscale_mat_uni <- function(x) {
y <- x
y[1, ] <- log(x[1, ])
y
}
#' Fitting Bivariate and univariate angular mixture models
#'
#' @importFrom gtools rdirichlet
#' @importFrom future.apply future_lapply
#'
#' @param model angular model whose mixtures are to be fitted. Available choices are \code{"vmsin", "vmcos"} and \code{"wnorm2"} for
#' bivariate data, and \code{"vm"} and \code{"wnorm"} for univariate data.
#' @param data data matrix (if bivarate, in which case it must have two columns) or vector. If outside, the values
#' are transformed into the scale \eqn{[0, 2\pi)}. *Note:* BAMBI cannot handle missing data. Missing values must
#' either be removed or properly imputed.
#' @param ncomp number of components in the mixture model. Must be a positive integer. vector values are not allowed.
#' If \code{comp == 1}, a single component model is fitted.
#' @param start_par list with elements \code{pmix} (ignored if \code{comp == 1}), together with \code{kappa1, kappa2, mu1} and \code{mu2},
#' for bivariate models, and \code{kappa} and \code{mu} for univariate models,
#' all being vectors of length same as \code{ncomp}.
#' These provides the starting values for the Markov chain; with \eqn{j}-th component of each vector corresponding to the \eqn{j}-th
#' component of the mixture distribution. If missing, the data is first clustered into \code{ncomp} groups either via k-means (after
#' projecting onto a unit sphere), or randomly, depending on \code{rand_start}, and then moment estimators for components are used as
#' the starting points. Note that a very wrong starting point can potentially lead the chain to get stuck at a wrong solution for thousands
#' of iterations. As such, we recommend using the default option, which is k-means followed by moment estimation.
#' @param cov.restrict Should there be any restriction on the covariance parameter for a bivariate model. Available choices are
#' \code{"POSITIVE", "NEGATIVE", "ZERO"} and "NONE". Note that \code{"ZERO"} fits a mixture with product components. Defaults to
#' \code{"NONE"}.
#' @param unimodal.component logical. Should each component in the mixture model be unimodal? Only used if \code{model} is either \code{"vmsin"}
#' or \code{"vmcos"}. Defaults to FALSE.
#' @param n.chains number of chains to run. Must be a positive integer.
#' @param chains_parallel logical. Should the chains be run in parallel? Defaluts to TRUE, and ignored if \code{n.chains} = 1.
#' Note that parallelization is implemented via \link{future_lapply} from package \code{future.apply} which
#' uses futures for this purpose, and thus provides a convenient way of parallelization across various OSs and computing environments.
#' However, a proper \link[future]{plan} must be set for the parallization before running the chain. Otherwise the chains will run sequentially.
#' @param method MCMC strategy to be used for the model paramters: \code{"hmc"} or \code{"rwmh"}.
#' @param perm_sampling logical. Should the permutation sampling algorithm of Fruhwirth-Schnatter (2001) be used?
#' If TRUE, at every iteration after burnin, once model parameters and mixing proportions are sampled,
#' a random permutation of 1, ..., ncomp is considered, and components are relabelled according
#' to this random permutation. This forced random label switchings may imporve the mixing rate of the chage. However, (automated) tuning
#' is very difficult with such a scheme, as there is no simple way of keeping track of the "original" component labels. This creates problem
#' with computing standard deviations of the generated model parameters, thus making the
#' scaling step used in tuning for \code{epsilon} or \code{paramscale} problematic as well. As such, \code{perm_sampling} is always turned
#' off during burn-in (even if \code{autotune = FALSE}), and turned on thereafter, if \code{TRUE}.
#' Defaults to and is set to \code{FALSE}.
#' @param int.displ absolute integer displacement for each coordinate for \code{wnorm} and \code{wnorm2} models (ignored otherwise). Default is 3.
#' Allowed minimum and maximum are 1 and 5 respectively.
#' @param epsilon,L tuning parameters for HMC; ignored if \code{method = "rwmh"}. \code{epsilon} (step-size) is a single number,
#' or a vector of size \code{2*ncomp} for univariate models and \code{5*ncomp} for bivariate models. Note that the "mass matrix"
#' in HMC is assumed to be identity. As such, \code{epsilon}'s corresponding to different model parameters need to be in proper scale for
#' optimal acceptance rate. Can be autotuned during burnin. See \code{autotune}.
#' \code{L} (leapfrog steps) is a positive integer or a vector of positive integers of length \code{n.chains}.
#' If multiple chains are used, we suggest same \code{L} values acorss different chains to make the chains as homogenous as possible.
#'
#' @param epsilon.random logical. Should \code{epsilon*delta}, where \code{delta} is a random
#' number between \code{(1-epsilon.incr, 1+epsilon.incr)} be used instead of \code{epsilon} at each iteration?
#' Ignored if \code{method = "rwmh"}.
#' @param L.random logical. Should a random integer between \code{L.orig/exp(L.incr)} and \code{L.orig*exp(L.incr)}be used instead as \code{L}
#' at each iteration? Ignored if \code{method = "rwmh"}. Defaults to \code{TRUE}.
#' @param L.incr amount of randomness incorporated in L if \code{L.random = TRUE}.
#' @param epsilon.incr amount of randomness incorporated in \code{epsilon} if \code{epsilon.random = TRUE}.
#' @param propscale tuning parameters for RWMH; a vector of size 5 (for bivariate models) or 2 (for univariate models) representing
#' the variances for the proposal normal densities
#' for the model parameters. Ignored if \code{method = "hmc"}. Can be autotuned during burnin. See \code{autotune}.
#' @param n.iter number of iterations for the Markov Chain.
## @param gam.loc,gam.scale location and scale (hyper-) parameters for the gamma prior for \code{kappa1} and \code{kappa2}. See
## \link{dgamma}. Defaults are \code{gam.loc = 0, gam.scale = 1000} that makes the prior non-informative.
#' @param pmix.alpha concentration parameter(s) for the Dirichlet prior for \code{pmix}. Must either be a positive real number, or a vector
#' with positive entries and of length \code{ncomp}. The default is \eqn{(r+r(r+1)/2)/2+3}, where \eqn{r} is 1 or 2 according as whether
#' the model is univariate or bivariate. Note that it is recommended to use larger \code{alpha} values to ensure the a good posterior behavior,
#' especially when \link{fit_incremental_angmix} is used for model selection, which handles overfitting in "let two component-specific parameters be
# identical", uses total number of components in the fitted model as the estimator for true component
#' size, and then penalizes for model complexity. See Fruhwirth-Schnatter (2011) for more details on this.
#' @param norm.var variance (hyper-) parameters in the normal prior for \code{log(kappa), log(kappa1), log(kappa2)} and \code{kappa3}. (Prior mean is zero).
#' Can be a vector. Default is 1000 that makes the prior non-informative.
#' @param burnin.prop proportion of iterations to used for burnin. Must be a be a number in [0, 1].
#' Default is 0.5.
#' @param thin thining size to be used. Must be a positive integer. If \code{thin = } n, then every nth iteration is reatained
#' in the final MCMC sample.
#' @param autotune logical. Should the Markov chain auto-tune the parameter \code{epsilon} (in HMC) or
#' \code{propscale} (in RWMH) during burn-in? Set to \code{TRUE} by default. An adaptive tuning strategy is implemented.
#' Here, at every 10th iteration during in burn-in, the acceptance ratio in the last \code{tune_ave_size}
#' iterations is calculated. Then the tuning parameter is decreased (increased) by a factor of
#' \code{1-tune.incr} (\code{1+tune.incr}) if the calculated acceptance rate
#' falls below (above) \code{accpt.prob.lower} (\code{accpt.prob.upper}). In addditon, when \code{iter} is a multiple of
#' \code{tune_ave_size}, \code{epsilon} for each model parameter is rescaled via the standard deviation of
#' the corresponding parameter over the past \code{tune_ave_size} iterations.
#' @param tune.prop proportion of *\code{burnin}* used to tune the parameters (\code{epsilon} in HMC and
#' \code{propscale} in RWMH). Must be a number between 0 and 1; defaults to 1. Ignored if \code{autotune == FALSE}.
#' @param show.progress logical. Should a progress bar be displayed?
#' @param accpt.prob.lower,accpt.prob.upper lower and upper limits of acceptance ratio to be maintained while tuning
#' during burn-in. Must be numbers between 0 and 1, which \code{accpt.prob.lower < accpt.prob.upper}. See \code{autotune}. Default to (0.6, 0,9) for HMC and (0.3, 0.5) for RWMH.
#' Ignored if \code{autotune = FALSE}.
#' @param tune.incr how much should the tuning parameter be increased or decreased at each step while tuning during burn-in?
#' Must be a number between 0 and 1. See \code{autotune}. Defaults to 0.05. Ignored if \code{autotune = FALSE}.
#' @param rand_start logical. Should a random starting clustering be used? Must be either a scalar, or a vector of length \code{ncomp},
#' one for each chain. Ignored if \code{start_par} is supplied. See \code{start_par} for more details. Defaults to \code{FALSE}.
#' @param tune_ave_size number previous iterations used to compute the acceptance rate while tuning in burn-in. Must be a positive
#' integer. Defaults to 100.
#' @param qrnd,n_qrnd Used only if \code{method="vmcos"}. See \link{dvmcos} for details.
#' @param kappa_upper,kappa_lower upper and lower bounds for the concentration and (absolute) association parameters. Must be a positive integers. Defaults to 150 and 1e-4,
#' and parameter with value above or below these limits rarely make sense in practice.
#' Warning: values much larger or smaller than the default are not recommended as they can cause numerical instability.
#' @param return_llik_contri logical. Should the log likelihood contribution of each data point for each MCMC iteration in each chain be returned? This makes
#' computation of \link{waic.angmcmc} and \link{loo.angmcmc} much faster. *Warning*: Depending on the length of data and \code{n.iter}, this can be
#' very memory intensive. We suggest setting \code{return_llik_contri = TRUE} only if \link{waic.angmcmc} and \link{loo.angmcmc} are aimed for. Defaults to
#' \code{FALSE}.
#' @param return_tune_param logical. Should the values of the tuning parameters used at each iteration in each chain be returned? Defaults to \code{FALSE}.
#' @param ... Unused.
#'
#' @note
#' Sampling is done in log scale for the concentration parameters (kappa, kappa1 and kappa2).
#'
#'
#' Parallelization is done by default when more than one chain is used,
#' but the chains can be run sequentially as well by setting
#' \code{chains_parallel = FALSE}. To retain reproducibility while running
#' multiple chains in parallel, the same RNG state is passed at the
#' beginning of each chain. This is done by specifying \code{future.seed = TRUE}
#' in \code{future.apply::future_lapply} call. Then at the beginning of the i-th
#' chain, before drawing any parameters, i-many Uniform(0, 1) random numbers are
#' generated using \code{runif(i)} (and then thrown away). This ensures that the
#' RNG states across chains prior to random generation of the parameters are
#' different, and hence, no two chains can become identical, even if they have
#' the same starting and tuning parameters. This, however creates a difference
#' between a \code{fit_angmix} call with multiple chains which is run sequentially
#' by setting \code{chains_parallel = FALSE}, and another which is run sequentially
#' because of a sequential \code{plan()} (or no \code{plan()}), with
#' \code{chains_parallel = TRUE}. In the former, different RNG states are passed at
#' the initiation of each chain.
#'
#'
#' @examples
#' # illustration only - more iterations needed for convergence
#' fit.vmsin.20 <- fit_angmix("vmsin", tim8,
#' ncomp = 3, n.iter = 20,
#' n.chains = 1
#' )
#' fit.vmsin.20
#'
#'
#' # Parallelization is implemented via future_lapply from the
#' # package future.apply. To parallelize, first provide a parallel
#' # plan(); otherwise the chains will run sequentially.
#' # Note that not all plan() might work on every OS, as they execute
#' # functions defined internally in fit_mixmodel. We suggest
#' # plan(multisession) which works on every OS.
#' \donttest{
#' library(future)
#' library(parallel)
#' # plan(multisession, gc = TRUE) # parallelize chains
#'
#' set.seed(1)
#' MC.fit <- fit_angmix("vmsin", tim8,
#' ncomp = 3, n.iter = 5000,
#' n.chains = 3
#' )
#'
#'
#' pointest(MC.fit)
#'
#' MC.fix <- fix_label(MC.fit)
#'
#' contour(MC.fit)
#' contour(MC.fix)
#' lpdtrace(MC.fit)
#' }
#'
#' @references
#' Fruhwirth-Schnatter, S. (2011). Label switching under model uncertainty. Mixtures: Estimation and Application, 213-239.
#'
#' Fruhwirth-Schnatter, S. (2001). Markov chain Monte Carlo estimation of classical and dynamic switching and mixture models. Journal of the American Statistical Association, 96(453), 194-209.
#'
#' @export
fit_angmix <- function(model = "vmsin",
data,
ncomp,
cov.restrict = "NONE",
unimodal.component = FALSE,
start_par = NULL,
rand_start = rep(FALSE, n.chains),
method = "hmc",
perm_sampling = FALSE,
n.chains = 3,
chains_parallel = TRUE,
return_llik_contri = FALSE,
int.displ = 3,
epsilon = 0.1,
L = 10,
epsilon.random = TRUE,
L.random = FALSE,
burnin.prop = 0.5,
tune.prop = 1,
thin = 1,
propscale = 0.05,
n.iter = 500,
# gam.loc = 0.001,
# gam.scale = 1000,
pmix.alpha = NULL,
norm.var = 1000,
autotune = TRUE,
show.progress = TRUE,
accpt.prob.upper,
accpt.prob.lower,
epsilon.incr = 0.05,
L.incr = 0.075,
tune.incr = 0.05,
tune_ave_size = 100,
kappa_upper = 150,
kappa_lower = 1e-4,
return_tune_param = FALSE,
qrnd = NULL,
n_qrnd = NULL, ...) {
# if(is.null(dim(data)) | !(mode(data) %in% c("list", "numeric", "data.frame") && ncol(data) == 2)) stop("non-compatible data")
dots <- list(...)
backward_compatible <- FALSE
backward_compatible <- dots$backward_compatible
if (is.null(backward_compatible)) {
backward_compatible <- FALSE
}
stopifnot(
is.logical(backward_compatible)
)
# progress.backend <- "tcltk"
# progressor <- function(...) NULL
progress.backend <- "tcltk"
signif_ <- if (backward_compatible) function(x, ...) x else function(x, ...) signif(x, 8)
if (length(model) > 1) stop("\'model\' must be a scalar")
if (model %in% c("vmsin", "vmcos", "wnorm2")) {
type <- "bi"
} else if (model %in% c("vm", "wnorm")) {
type <- "uni"
} else {
stop("non-compatible model")
}
if (missing(ncomp)) {
stop("\'ncomp\' is missing, with no default.")
}
if (length(ncomp) > 1) {
stop("\'ncomp\' must be a scalar")
}
if (any(length(n.chains) > 1, length(n.chains) < 1, n.chains == 0)) {
stop("Invalid n.chains")
}
zero_cov <- FALSE
if (type == "bi") {
if (!(is.matrix(data) | is.data.frame(data))) {
stop("\'data\' must be a two column matrix for model = \'vmsin\', \'vmcos\' and \'wnorm2\'")
}
if (ncol(data) != 2) {
stop("\'data\' must be a two column matrix for model = \'vmsin\', \'vmcos\' and \'wnorm2\'")
}
if (length(cov.restrict) > 1 |
!cov.restrict %in% c("NONE", "POSITIVE", "NEGATIVE", "ZERO")) {
stop("cov.restrict must be one of \"NONE\", \"POSITIVE\", \"NEGATIVE\" and \"ZERO\"")
}
if (any(is.na(data))) {
stop(paste(
"\'data\' contains missing value(s). BAMBI cannot handle missing data.\n",
"Either remove them or properly impute them."
))
}
data.rad <- rm_NA_rad(data)
n.data <- nrow(data.rad)
npar_1_comp <- 5
par.names <- c("kappa1", "kappa2", "kappa3", "mu1", "mu2")
par_lower <- replicate(ncomp, c(kappa_lower, kappa_lower, -kappa_upper, 0, 0))
par_upper <- replicate(ncomp, c(
kappa_upper, kappa_upper,
kappa_upper, 2 * pi, 2 * pi
))
if (cov.restrict == "POSITIVE") {
par_lower[3, ] <- 0
} else if (cov.restrict == "NEGATIVE") {
par_upper[3, ] <- 0
} else if (cov.restrict == "ZERO") {
par_lower[3, ] <- par_upper[3, ] <- 0
zero_cov <- TRUE
}
if (missing(pmix.alpha)) {
pmix.alpha <- 5.5
}
} else {
if ((is.matrix(data) | is.data.frame(data))) {
data <- as.vector(as.matrix(data))
}
if (!is.numeric(data)) {
stop("\'data\' must be a vector for \'model\' = \'vm\' and \'wnorm\'")
}
data <- as.numeric(data)
data.rad <- rm_NA_rad(data)
n.data <- length(data.rad)
npar_1_comp <- 2
par.names <- c("kappa", "mu")
par_lower <- replicate(ncomp, c(kappa_lower, 0))
par_upper <- replicate(ncomp, c(kappa_upper, 2 * pi))
if (missing(pmix.alpha)) {
pmix.alpha <- 4
}
}
if (any(
length(perm_sampling) > 1, length(perm_sampling) < 1,
!is.logical(perm_sampling)
)) {
stop("Invalid perm_sampling")
}
if (any(
length(chains_parallel) > 1, length(chains_parallel) < 1,
!is.logical(chains_parallel)
)) {
stop("Invalid chains_parallel")
}
if (any(
length(chains_parallel) > 1, length(chains_parallel) < 1,
!is.logical(chains_parallel)
)) {
stop("Invalid chains_parallel")
}
if (any(
length(autotune) > 1, length(autotune) < 1,
!is.logical(autotune)
)) {
stop("Invalid autotune")
}
if (any(
length(unimodal.component) > 1, length(unimodal.component) < 1,
!is.logical(unimodal.component)
)) {
stop("Invalid unimodal.component")
}
if (any(
length(return_llik_contri) > 1, length(return_llik_contri) < 1,
!is.logical(return_llik_contri)
)) {
stop("Invalid return_llik_contri")
}
if (any(
length(return_tune_param) > 1, length(return_tune_param) < 1,
!is.logical(return_tune_param)
)) {
stop("Invalid return_tune_param")
}
if (length(pmix.alpha) == 1) {
pmix.alpha <- rep(pmix.alpha, ncomp)
} else if (length(pmix.alpha) != ncomp) {
stop("length(pmix.alpha) and ncomp differ")
}
if (n.chains < 1) {
stop("\'n.chains\' must be a positive integer")
}
if (n.chains == 1) {
chains_parallel <- FALSE
}
if (length(burnin.prop) != 1 || burnin.prop < 0 || burnin.prop >= 1) {
stop("\"burnin.prop\" must be a number in [0, 1)")
}
if (thin < 1) {
stop("\"thin\" must be a positive integer")
}
if (tune.incr <= 0 | tune.incr >= 1) {
stop("\'tune.incr\' must be between 0 and 1")
}
n.burnin <- ceiling(burnin.prop * n.iter)
if (length(tune.prop) != 1 || tune.prop < 0 || tune.prop > 1) {
stop("\"tune.prop\" must be in [0, 1]")
}
iter.tune <- ceiling(burnin.prop * tune.prop * n.iter)
if (iter.tune > n.burnin) {
iter.tune <- n.burnin
}
n.iter.final <- n.iter - n.burnin
burnin_iter <- seq_len(n.burnin)
thin <- round(thin)
thin_filter <- c(TRUE, rep(FALSE, thin - 1))
final_iter_set <- (seq_len(n.iter))[-burnin_iter][thin_filter]
# gam.rate <- 1/gam.scale
curr.model <- model
if (length(rand_start) == 1) {
rand_start <- rep(rand_start, n.chains)
}
if (method == "hmc") {
if (any(L < 1)) {
stop("\'L\' must be a positive integer")
}
L <- ceiling(L)
if (length(L) == 1) {
L <- rep(L, n.chains)
} else if (length(L) != n.chains) {
stop("\'L\' must be a vector of length \'n.chains\'")
} else {
rand_start <- rep(rand_start, n.chains)
}
if (missing(accpt.prob.upper)) {
accpt.prob.upper <- 0.9
}
if (missing(accpt.prob.lower)) {
accpt.prob.lower <- if (model %in% c("wnorm", "wnorm2")) 0.58 else 0.6
}
if (length(epsilon) == ncomp * npar_1_comp) {
tune_param <- matrix(c(epsilon), npar_1_comp, ncomp)
} else if (length(epsilon) == 1) {
tune_param <- matrix(epsilon, npar_1_comp, ncomp)
} else {
stop("epsilon must either be a scalar or a vector of length 2*ncomp (univariate) or 5*ncomp (bivariate)")
}
}
#
# if (method == "hmc" & sum(rand_start) == 0 & n.chains > 1 & chains_parallel) {
# if (length(L) == 1 | max(L) == min(L)) {
# warning(paste("same L accross multiple chains while running them in",
# "parllel with rand_start = FALSE will just make",
# "identical copies of the same chain. L is changed"))
# L <- seq(ceiling(L/2), ceiling(2*L), length.out = n.chains)
# }
# }
if (grepl(method, "rwmh")) # using rwmh
{
if (missing(accpt.prob.upper)) {
accpt.prob.upper <- 0.5
}
if (missing(accpt.prob.lower)) {
accpt.prob.lower <- 0.3
}
if (length(propscale) == ncomp * npar_1_comp) {
tune_param <- matrix(c(propscale), npar_1_comp, ncomp)
} else if (length(propscale) == 1) {
tune_param <- matrix(propscale, npar_1_comp, ncomp)
} else {
stop("propscale must either be a scalar or a vector of length 2*ncomp (univariate) or 5*ncomp (bivariate)")
}
}
nterms <- tune_ave_size
if (!model %in% c("wnorm", "wnorm2")) {
int.displ <- omega.2pi <- NULL
}
if (model != "vmcos") {
qrnd_grid <- NULL
n_qrnd <- NULL
}
# Now rename the model specific compiled llik and grad functions
if (model == "vmsin") {
if (unimodal.component) {
dep_cov <- TRUE
dep_cov_type <- "vmsin_unimodal"
} else {
dep_cov <- FALSE
dep_cov_type <- NULL
}
# # in log scale
# lpd_grad_model_indep_1comp <- function(data, par_vec_lscale,
# obs_group, n.clus) {
# par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
# lpd_grad <- matrix(NA, 6, 1)
# if (n.clus > 0) {
# lpd_grad <- (grad_llik_vmsin_C(data[obs_group, , drop=FALSE],
# par_vec) +
# c( # grad for lprior
# (gam.loc - 1)/par_vec[1:2] - gam.rate, -par_vec[3]/norm.var, 0, 0,
# # lprior
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# )) * c(par_vec[1:2], rep(1, 4))
# } else {
# lpd_grad <-
# c( # grad for lprior
# (gam.loc - 1)/par_vec[1:2] - gam.rate, -par_vec[3]/norm.var, 0, 0,
# # lprior
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# ) *
# c(par_vec[1:2], rep(1, 4))
# }
#
# list(lpr = (lpd_grad[6]), grad = lpd_grad[1:5])
# }
#
#
# lpd_model_indep_1comp <- function(data, par_vec_lscale, obs_group,
# n.clus) {
# par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
# if (n.clus > 0) {
# # llik + prior
# res <-
# llik_vmsin_one_comp(data[obs_group, , drop=FALSE], par_vec,
# log(const_vmsin(par_vec[1],
# par_vec[2], par_vec[3]))) +
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# } else{
# # only prior
# res <-
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# }
# res
# }
# in log scale
lpd_grad_model_indep_1comp <- function(data, par_vec_lscale,
obs_group, n.clus) {
par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
lpd_grad <- matrix(NA, 6, 1)
if (n.clus > 0) {
lpd_grad <- signif_(suppressWarnings(
grad_llik_vmsin_C(
data[obs_group, , drop = FALSE],
par_vec
)
) *
c(par_vec[1:2], rep(1, 4)) +
c( # grad for lprior
-par_vec_lscale[1:3] / norm.var, 0, 0,
# lprior
-0.5 * sum(par_vec_lscale[1:3]^2 / norm.var)
))
} else {
lpd_grad[] <-
signif_(c( # grad for lprior
-par_vec_lscale[1:3] / norm.var, 0, 0,
# lprior
-0.5 * sum(par_vec_lscale[1:3]^2 / norm.var)
))
}
list(lpr = (lpd_grad[6]), grad = lpd_grad[1:5])
}
lpd_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
if (n.clus > 0) {
# llik + prior
res <-
signif_(llik_vmsin_one_comp(
data[obs_group, , drop = FALSE], par_vec,
log(const_vmsin(
par_vec[1],
par_vec[2], par_vec[3]
))
) +
0.5 * sum(-par_vec_lscale[1:3]^2 / norm.var))
} else {
# only prior
res <-
signif_(0.5 * sum(-par_vec_lscale[1:3]^2 / norm.var))
}
res
}
llik_model_contri <- function(data, par_mat, pi_mix) {
signif_(
llik_vmsin_contri_C(
data, par_mat, pi_mix,
signif_(log_const_vmsin_all(par_mat))
)
)
}
mem_p_model <- function(data, par_mat, pi_mix) {
signif_(
mem_p_sin(
data, par_mat, pi_mix,
signif_(log_const_vmsin_all(par_mat)), 1
)
)
}
} else if (model == "vmcos") {
ell <- list(qrnd_grid = qrnd, n_qrnd = n_qrnd)
if (!is.null(ell$qrnd)) {
qrnd_grid <- ell$qrnd
dim_qrnd <- dim(qrnd_grid)
if (!is.matrix(qrnd_grid) | is.null(dim_qrnd) |
dim_qrnd[2] != 2) {
stop("qrnd_grid must be a two column matrix")
}
n_qrnd <- dim_qrnd[1]
} else if (!is.null(ell$n_qrnd)) {
n_qrnd <- round(ell$n_qrnd)
if (n_qrnd < 1) {
stop("n_qrnd must be a positive integer")
}
qrnd_grid <- sobol(n_qrnd, 2, FALSE)
} else {
n_qrnd <- 1e4
qrnd_grid <- sobol(n_qrnd, 2, FALSE)
}
if (unimodal.component) {
dep_cov <- TRUE
dep_cov_type <- "vmcos_unimodal"
} else {
dep_cov <- FALSE
dep_cov_type <- NULL
}
# # in log scale
# lpd_grad_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
# par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
# lpd_grad <- matrix(NA, 6, 1)
# if (n.clus > 0) {
# lpd_grad[] <- (grad_llik_vmcos_C(data[obs_group, , drop=FALSE],
# par_vec[], qrnd_grid) +
# c( # grad for lprior
# (gam.loc - 1)/par_vec[1:2] - gam.rate, -par_vec[3]/norm.var, 0, 0,
# # lprior
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var)
# ) *
# c(par_vec[1:2], rep(1, 4))
# } else {
# lpd_grad[] <-
# c( # grad for lprior
# (gam.loc - 1)/par_vec[1:2] - gam.rate, -par_vec[3]/norm.var, 0, 0,
# # lprior
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# ) *
# c(par_vec[1:2], rep(1, 4))
# }
# list(lpr = (lpd_grad[6]), grad = lpd_grad[1:5])
# }
#
# lpd_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
#
# par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
# if (n.clus > 0) {
# # llik + prior
# res <-
# llik_vmcos_one_comp(data[obs_group, , drop=FALSE], par_vec[],
# log(const_vmcos(par_vec[1],
# par_vec[2],
# par_vec[3],
# qrnd_grid))) +
# sum((gam.loc - 1)*log(par_vec[1:2])-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# } else{
# # only prior
# res <-
# sum((gam.loc - 1)*log(par_vec[1:2])-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# }
# res
# }
# in log scale
lpd_grad_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
lpd_grad <- matrix(NA, 6, 1)
if (n.clus > 0) {
lpd_grad[] <- signif_(suppressWarnings(
grad_llik_vmcos_C(
data[obs_group, , drop = FALSE],
par_vec[], qrnd_grid
)
) *
c(par_vec[1:2], rep(1, 4)) +
c( # grad for lprior
-par_vec_lscale[1:3] / norm.var, 0, 0,
# lprior
-0.5 * sum(par_vec_lscale[1:3]^2 / norm.var)
))
} else {
lpd_grad[] <-
signif_(c( # grad for lprior
-par_vec_lscale[1:3] / norm.var, 0, 0,
# lprior
-0.5 * sum(par_vec_lscale[1:3]^2 / norm.var)
))
}
list(lpr = (lpd_grad[6]), grad = lpd_grad[1:5])
}
lpd_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
if (n.clus > 0) {
# llik + prior
res <-
signif_(llik_vmcos_one_comp(
data[obs_group, , drop = FALSE], par_vec[],
log(const_vmcos(
par_vec[1],
par_vec[2],
par_vec[3],
qrnd_grid
))
) -
0.5 * sum(par_vec_lscale[1:3]^2 / norm.var))
} else {
# only prior
res <-
signif_(-0.5 * sum(par_vec_lscale[1:3]^2 / norm.var))
}
res
}
llik_model_contri <- function(data, par_mat, pi_mix) {
signif_(
llik_vmcos_contri_C(
data, par_mat, pi_mix,
signif_(log_const_vmcos_all(par_mat, qrnd_grid))
)
)
}
mem_p_model <- function(data, par_mat, pi_mix) {
signif_(
mem_p_cos(
data, par_mat, pi_mix,
signif_(log_const_vmcos_all(par_mat, qrnd_grid))
)
)
}
} else if (model == "wnorm2") {
dep_cov <- TRUE
dep_cov_type <- "wnorm2_bound"
if (int.displ >= 5) {
int.displ <- 5
} else if (int.displ <= 1) int.displ <- 1
int.displ <- floor(int.displ)
omega.2pi.all <- expand.grid(-int.displ:int.displ, -int.displ:int.displ) * (2 * pi) # 2pi * integer displacements
omega.2pi <- as.matrix(omega.2pi.all)
# lpd_grad_model_indep <- function(data, par_mat, obs_group, n.clus) {
# lpd_grad <- matrix(NA, 6, ncomp)
# for(j in 1:ncomp) {
# if (n.clus[j] > 0) {
# lpd_grad[, j] <- grad_llik_wnorm2_C(data[obs_group[[j]], , drop=FALSE],
# par_mat[, j], omega.2pi) +
# c( # grad for lprior
# (gam.loc - 1)/par_mat[1:2, j] - gam.rate, -par_mat[3, j]/norm.var, 0, 0,
# # lprior
# sum((gam.loc - 1)*log(par_mat[1:2, j])-
# gam.rate*par_mat[1:2, j]) - 0.5*par_mat[3, j]^2/norm.var
# )
# } else {
# lpd_grad[, j] <-
# c( # grad for lprior
# (gam.loc - 1)/par_mat[1:2, j] - gam.rate, -par_mat[3, j]/norm.var, 0, 0,
# # lprior
# sum((gam.loc - 1)*log(par_mat[1:2, j])-
# gam.rate*par_mat[1:2, j]) - 0.5*par_mat[3, j]^2/norm.var
# )
# }
# }
# list(lpr = sum(lpd_grad[6, ]), grad = lpd_grad[1:5, ])
# }
# lpd_model_indep <- function(data, par_mat, obs_group, n.clus) {
# res <- 0
# for(j in 1:ncomp) {
# if (n.clus[j] > 0) {
# # llik + prior
# res <- res +
# llik_wnorm2_one_comp(data[obs_group[[j]], , drop=FALSE], par_mat[, j],
# l_const_wnorm2(par_mat[, j]),
# omega.2pi) +
# sum((gam.loc - 1)*log(par_mat[1:2, j])-
# gam.rate*par_mat[1:2, j]) - 0.5*par_mat[3, j]^2/norm.var
# } else{
# # only prior
# res <- res +
# sum((gam.loc - 1)*log(par_mat[1:2, j])-
# gam.rate*par_mat[1:2, j]) - 0.5*par_mat[3, j]^2/norm.var
# }
# }
# unname(res)
# }
# # in log scale
# lpd_grad_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
# par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
# lpd_grad <- matrix(NA, 6, 1)
# if (n.clus > 0) {
# lpd_grad[] <- (grad_llik_wnorm2_C(data[obs_group, , drop=FALSE],
# par_vec[], omega.2pi) +
# c( # grad for lprior
# (gam.loc - 1)/par_vec[1:2] - gam.rate, -par_vec[3]/norm.var, 0, 0,
# # lprior
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# )) *
# c(par_vec[1:2], rep(1, 4))
# } else {
# lpd_grad[] <-
# c( # grad for lprior
# (gam.loc - 1)/par_vec[1:2] - gam.rate, -par_vec[3]/norm.var, 0, 0,
# # lprior
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# ) *
# c(par_vec[1:2], rep(1, 4))
# }
# list(lpr = sum(lpd_grad[6]), grad = lpd_grad[1:5])
# }
#
# lpd_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
# par_vec <- c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5])
# if (n.clus > 0) {
# # llik + prior
# res <-
# llik_wnorm2_one_comp(data[obs_group, , drop=FALSE], par_vec[],
# l_const_wnorm2(par_vec[]),
# omega.2pi) +
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# } else{
# # only prior
# res <-
# sum((gam.loc - 1)*par_vec_lscale[1:2]-
# gam.rate*par_vec[1:2]) - 0.5*par_vec[3]^2/norm.var
# }
# unname(res)
# }
# in log scale
lpd_grad_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
par_vec <- signif_(c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5]))
lpd_grad <- matrix(NA, 6, 1)
if (n.clus > 0) {
lpd_grad[] <- signif_(
grad_llik_wnorm2_C(
data[obs_group, , drop = FALSE],
par_vec[], omega.2pi
) *
c(par_vec[1:2], rep(1, 4)) +
c( # grad for lprior
-par_vec_lscale[1:3] / norm.var, 0, 0,
# lprior
-0.5 * sum(par_vec_lscale[1:3]^2 / norm.var)
)
)
} else {
lpd_grad[] <-
signif_(c( # grad for lprior
-par_vec_lscale[1:3] / norm.var, 0, 0,
# lprior
-0.5 * sum(par_vec_lscale[1:3]^2 / norm.var)
))
}
list(lpr = (lpd_grad[6]), grad = lpd_grad[1:5])
}
lpd_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
par_vec <- signif_(c(exp(par_vec_lscale[1:2]), par_vec_lscale[3:5]))
if (n.clus > 0) {
# llik + prior
res <-
signif_(llik_wnorm2_one_comp(
data[obs_group, , drop = FALSE], par_vec[],
l_const_wnorm2(par_vec[]),
omega.2pi
) -
0.5 * sum(par_vec_lscale[1:3]^2 / norm.var))
} else {
# only prior
res <-
signif_(-0.5 * sum(par_vec_lscale[1:3]^2 / norm.var))
}
res
}
llik_model_contri <- function(data, par_mat, pi_mix) {
signif_(
llik_wnorm2_contri_C(
data, par_mat, pi_mix,
signif_(log_const_wnorm2_all(par_mat)),
omega.2pi
)
)
}
mem_p_model <- function(data, par_mat, pi_mix) {
signif_(
mem_p_wnorm2(
data, par_mat, pi_mix,
signif_(log_const_wnorm2_all(par_mat)),
omega.2pi
)
)
}
} else if (model == "vm") {
dep_cov <- FALSE
# lpd_grad_model_indep <- function(data, par_mat, obs_group, n.clus) {
# lpd_grad <- matrix(NA, 3, ncomp)
# for(j in 1:ncomp) {
# if (n.clus[j] > 0) {
# lpd_grad[, j] <- grad_llik_univm_C(data[obs_group[[j]]],
# par_mat[, j]) +
# c( # grad for lprior
# (gam.loc - 1)/par_mat[1, j] - gam.rate, 0,
# # lprior
# (gam.loc - 1)*log(par_mat[1, j])-
# gam.rate*par_mat[1, j]
# )
#
# } else {
# lpd_grad[, j] <-
# c( # grad for lprior
# (gam.loc - 1)/par_mat[1, j] - gam.rate, 0,
# # lprior
# (gam.loc - 1)*log(par_mat[1, j])-
# gam.rate*par_mat[1, j]
# )
# }
# }
# list(lpr = sum(lpd_grad[3, ]), grad = lpd_grad[1:2, ])
# }
#
# lpd_model_indep <- function(data, par_mat, obs_group, n.clus) {
# res <- 0
# for(j in 1:ncomp) {
# if (n.clus[j] > 0) {
# # llik + prior
# res <- res +
# llik_univm_one_comp(data[obs_group[[j]]], par_mat[, j],
# log(const_univm(par_mat[1, j]))) +
# (gam.loc - 1)*log(par_mat[1, j]) - gam.rate*par_mat[1, j]
#
# } else{
# # only prior
# res <- res +
# (gam.loc - 1)*log(par_mat[1, j]) - gam.rate*par_mat[1, j]
# }
# }
# unname(res)
# }
# lpd_grad_model_indep_1comp <- function(data, par_vec, obs_group, n.clus) {
# lpd_grad <- matrix(NA, 3, 1)
# if (n.clus > 0) {
# lpd_grad[] <- grad_llik_univm_C(data[obs_group],
# par_vec[]) +
# c( # grad for lprior
# (gam.loc - 1)/par_vec[1] - gam.rate, 0,
# # lprior
# (gam.loc - 1)*log(par_vec[1])-
# gam.rate*par_vec[1]
# )
#
# } else {
# lpd_grad[] <-
# c( # grad for lprior
# (gam.loc - 1)/par_vec[1] - gam.rate, 0,
# # lprior
# (gam.loc - 1)*log(par_vec[1])-
# gam.rate*par_vec[1]
# )
# }
# list(lpr = (lpd_grad[3 ]), grad = lpd_grad[1:2 ])
# }
#
# in log scale
lpd_grad_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
lpd_grad <- matrix(NA, 3, 1)
par_vec <- c(exp(par_vec_lscale[1]), par_vec_lscale[2])
if (n.clus > 0) {
lpd_grad[] <- signif_(suppressWarnings(
grad_llik_univm_C(
data[obs_group],
par_vec[]
)
) * c(par_vec[1], 1, 1) +
c( # grad for lprior
-par_vec_lscale[1] / norm.var, 0,
# lprior
-0.5 * sum(par_vec_lscale[1]^2 / norm.var)
))
} else {
lpd_grad[] <-
signif_(c( # grad for lprior
-par_vec_lscale[1] / norm.var, 0,
# lprior
-0.5 * sum(par_vec_lscale[1]^2 / norm.var)
))
}
list(lpr = (lpd_grad[3]), grad = lpd_grad[1:2])
}
lpd_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
par_vec <- c(exp(par_vec_lscale[1]), par_vec_lscale[2])
if (n.clus > 0) {
# llik + prior
res <-
signif_(
llik_univm_one_comp(
data[obs_group], par_vec[],
log(const_univm(par_vec[1]))
) -
0.5 * sum(par_vec_lscale[1]^2 / norm.var)
)
} else {
# only prior
res <-
signif_(-0.5 * sum(par_vec_lscale[1]^2 / norm.var))
}
unname(res)
}
llik_model_contri <- function(data, par_mat, pi_mix) {
signif_(
llik_univm_contri_C(
data, par_mat, pi_mix,
signif_(log_const_univm_all(par_mat))
)
)
}
mem_p_model <- function(data, par_mat, pi_mix) {
signif_(
mem_p_univm(
data, par_mat, pi_mix,
signif_(log_const_univm_all(par_mat))
)
)
}
}
# else if (model == "wnorm")
else {
dep_cov <- FALSE
if (int.displ >= 5) {
int.displ <- 5
} else if (int.displ <= 1) int.displ <- 1
int.displ <- floor(int.displ)
omega.2pi <- (-int.displ):int.displ * (2 * pi) # 2pi * 1d integer displacements
#
# lpd_grad_model_indep <- function(data, par_mat, obs_group, n.clus) {
# lpd_grad <- matrix(NA, 3, ncomp)
# for(j in 1:ncomp) {
# if (n.clus[j] > 0) {
# lpd_grad[, j] <- grad_llik_uniwnorm_C(data[obs_group[[j]]],
# par_mat[, j], omega.2pi) +
# c( # grad for lprior
# (gam.loc - 1)/par_mat[1, j] - gam.rate, 0,
# # lprior
# (gam.loc - 1)*log(par_mat[1, j])-
# gam.rate*par_mat[1, j]
# )
#
# } else {
# lpd_grad[, j] <-
# c( # grad for lprior
# (gam.loc - 1)/par_mat[1, j] - gam.rate, 0,
# # lprior
# (gam.loc - 1)*log(par_mat[1, j])-
# gam.rate*par_mat[1, j]
# )
# }
# }
# list(lpr = sum(lpd_grad[3, ]), grad = lpd_grad[1:2, ])
# }
#
# lpd_model_indep <- function(data, par_mat, obs_group, n.clus) {
# res <- 0
# for(j in 1:ncomp) {
# if (n.clus[j] > 0) {
# # llik + prior
# res <- res +
# llik_uniwnorm_one_comp(data[obs_group[[j]]], par_mat[, j],
# l_const_uniwnorm(par_mat[1, j]),
# omega.2pi) +
# (gam.loc - 1)*log(par_mat[1, j]) - gam.rate*par_mat[1, j]
#
# } else{
# # only prior
# res <- res +
# (gam.loc - 1)*log(par_mat[1, j]) - gam.rate*par_mat[1, j]
# }
# }
# unname(res)
# }
lpd_grad_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
lpd_grad <- matrix(NA, 3, 1)
par_vec <- c(exp(par_vec_lscale[1]), par_vec_lscale[2])
if (n.clus > 0) {
lpd_grad[] <- signif_(grad_llik_uniwnorm_C(
data[obs_group],
par_vec[], omega.2pi
) * c(par_vec[1], 1, 1) +
c( # grad for lprior
-par_vec_lscale[1] / norm.var, 0,
# lprior
-0.5 * sum(par_vec_lscale[1]^2 / norm.var)
))
} else {
lpd_grad[] <-
signif_(c( # grad for lprior
-par_vec_lscale[1] / norm.var, 0,
# lprior
-0.5 * sum(par_vec_lscale[1]^2 / norm.var)
))
}
list(lpr = (lpd_grad[3]), grad = lpd_grad[1:2])
}
lpd_model_indep_1comp <- function(data, par_vec_lscale, obs_group, n.clus) {
par_vec <- c(exp(par_vec_lscale[1]), par_vec_lscale[2])
if (n.clus > 0) {
# llik + prior
res <-
signif_(llik_uniwnorm_one_comp(
data[obs_group], par_vec[],
l_const_uniwnorm(par_vec[1]),
omega.2pi
) -
0.5 * sum(par_vec_lscale[1]^2 / norm.var))
} else {
# only prior
res <-
signif_(-0.5 * sum(par_vec_lscale[1]^2 / norm.var))
}
unname(res)
}
llik_model_contri <- function(data, par_mat, pi_mix) {
signif_(
llik_uniwnorm_contri_C(
data, par_mat, pi_mix,
signif_(log_const_uniwnorm_all(par_mat)),
omega.2pi
)
)
}
mem_p_model <- function(data, par_mat, pi_mix) {
signif_(
mem_p_uniwnorm(
data, par_mat, pi_mix,
signif_(log_const_uniwnorm_all(par_mat)),
omega.2pi
)
)
}
# grad_llik_uniwnorm_R <- function(data, par) {
# gr <- numDeriv::grad(function(par) llik_uniwnorm_one_comp(data, par,
# l_const_uniwnorm(par[1]),
# omega.2pi),
# par)
# gr[3] <- llik_uniwnorm_one_comp(data, par,
# l_const_uniwnorm(par[1]), omega.2pi)
# gr
# }
# llik_uniwnorm_full(data.rad, par.mat, pi.mix, log_c = log_const_uniwnorm_all(par.mat), omega.2pi)
#
#
# grad_llik_uniwnorm_C(data.rad, par.mat[, 1], omega.2pi)
# grad_llik_uniwnorm_R(data.rad, par.mat[, 1])
}
if (!is.null(start_par)) {
if (!is.list(start_par)) {
stop("start_par must be a list")
}
if (!is.list(start_par[[1]])) {
if (length(setdiff(c(par.names, "pmix"), names(start_par))) > 0) {
stop("start_par does not have all parameters")
}
start_par <- lapply(1:n.chains, function(ii) start_par)
} else {
if (length(start_par) != n.chains) {
stop("length(start_par) must be either 1 or equal to n.chains")
}
for (jj in 1:n.chains) {
if (length(setdiff(c(par.names, "pmix"), names(start_par[[jj]]))) > 0) {
stop(paste0("start_par[[", jj, "]] does not have all parameters"))
}
}
}
}
if (any(is.null(start_par), !is.list(start_par[[1]])) &
all(!rand_start)) {
starting_1 <- process_startpar(
start_par,
data.rad,
ncomp,
model,
FALSE
)
starting <- lapply(1:n.chains, function(j) starting_1)
} else {
starting <- lapply(
1:n.chains,
function(j) {
process_startpar(
start_par[[j]],
data.rad,
ncomp,
model,
rand_start[j]
)
}
)
}
run_MC <- function(starting, L, chain_no) {
# just to change the RNG state across chains, so that
# no two c chains turn out to be identical
change_rng_state <- runif(chain_no)
if (method == "hmc") # using hmc
{
L_vec <- rep(L, n.iter)
if (L.random) {
L.orig <- L
}
}
starting$par.mat[abs(starting$par.mat) >= kappa_upper / 2] <- kappa_upper / 2
starting$par.mat[abs(starting$par.mat) <= 2 * kappa_lower] <- 2 * kappa_lower
par.mat.all <- array(0, dim = c(npar_1_comp, ncomp, n.iter))
par.mat_lscale.all <- par.mat.all
pi.mix.all <- matrix(1, nrow = ncomp, ncol = n.iter)
llik.all <- lprior.all <- lpd.all <- rep(-Inf, (n.iter))
accpt.par.mat.all <- matrix(NA, ncomp, n.iter)
modelpar.names <- par.names
clus_ind.all <- matrix(1, nrow = n.data, ncol = n.iter)
if (return_llik_contri) {
llik_contri_all <- matrix(1, nrow = n.data, ncol = n.iter)
} else {
llik_contri_all <- NULL
}
mem_prob_all <- array(1, dim = c(n.data, ncomp, n.iter))
epsilon_ave <- NULL
L_ave <- NULL
propscale_final <- NULL
pi.mix <- starting$pi.mix
par.mat <- as.matrix(starting$par.mat)
if (zero_cov) {
par.mat[3, ] <- 0
} else if (cov.restrict == "POSITIVE") {
par.mat[3, ] <- pmax(par.mat[3, ], 0)
} else if (cov.restrict == "NEGATIVE") {
par.mat[3, ] <- pmin(par.mat[3, ], 0)
}
if (type == "bi") {
par.mat_lscale <- find_lscale_mat(par.mat)
par_lower_lscale <- find_lscale_mat(par_lower)
par_upper_lscale <- find_lscale_mat(par_upper)
} else {
par.mat_lscale <- find_lscale_mat_uni(par.mat)
par_lower_lscale <- find_lscale_mat_uni(par_lower)
par_upper_lscale <- find_lscale_mat_uni(par_upper)
}
# browser()
if (ncomp == 1) {
perm_sampling <- FALSE
clus.ind <- clus.ind_curr <- rep(1, n.data)
obs_group <- list(1:n.data)
n.clus <- n.data
post.wt <- matrix(1, n.data, 1)
pi.mix <- 1
}
tcltk_fail <- FALSE
if (show.progress & !exists("pb") & progress.backend == "tcltk") {
pb <- tryCatch(tcltk::tkProgressBar(
title = paste("Chain", chain_no),
label = "Initializing...",
min = 1, max = n.iter
),
error = function(e) e
)
if (is(pb, "error") | is(pb, "warning")) {
show.progress <- FALSE
tcltk_fail <- TRUE
}
} else if (show.progress & !exists("pb") & progress.backend == "txt") {
pb <- utils::txtProgressBar(title = paste("Chain", chain_no),
label = "Initializing...",
min = 1, max = n.iter, style = 3)
tcltk_fail <- FALSE
}
if (tcltk_fail) {
msg <- paste(
"{tcltk} could not be loaded; 'show.progress' was set to FALSE."#,
# "Consider setting progress.backend = \"txt\" for a text-based",
# "progress bar. See ?fit_angmix for more details."
)
message(msg)
show.progress <- FALSE
}
par.mat.all.order <- par.mat.all
ntunes_up <- ntunes_down <- rep(0, ncomp)
# tune_iter_no <- c()
# ave_accpt_all <- c()
tune_param_all <- matrix(NA, length(tune_param), n.iter)
tune_status <- matrix(0, ncomp, n.iter)
# Run the Markov chain
for (iter in seq_len(n.iter)) {
# browser()
if (method == "hmc") {
if (epsilon.random) {
tune_param_final <- tune_param * runif(1, 1 - epsilon.incr, 1 + epsilon.incr)
} else {
tune_param_final <- tune_param
}
if (L.random) {
L_vec[iter] <- L <- sample(ceiling(L.orig / exp(L.incr)):ceiling(L.orig * exp(L.incr)), 1)
}
}
#----------------------------------------------------------------------------------
# generating mixture proportions if ncomp > 1
#----------------------------------------------------------------------------------
if (ncomp > 1) {
post.wt <- mem_p_model(data.rad, par.mat, pi.mix)
clus.ind <- cID(post.wt, ncomp, runif(n.data))
# clus.ind <- apply(post.wt, 1, function(x) which.max(rmultinom(n = 1, size = 1, prob = x)))
# n.clus <- tabulate(clus.ind_curr, nbins = ncomp)
obs_group <- lapply(1:ncomp, function(j) which(clus.ind == j))
n.clus <- listLen(obs_group)
pi.mix <-
as.numeric(rdirichlet(1, (pmix.alpha + n.clus))) # new mixture proportions
}
#----------------------------------------------------------------------------------
# generating par.mat
#----------------------------------------------------------------------------------
if (type == "bi") {
if (method == "hmc") {
# hmc_curr <-
# bounded_hmc_biv(lpr_grad =
# function(par_mat)
# lpd_grad_model_indep(data.rad, par_mat,
# obs_group, n.clus),
# init = par.mat,
# lower = par_lower,
# upper = par_upper,
# dep_cov = dep_cov,
# dep_cov_type = dep_cov_type,
# zero_cov = zero_cov,
# nsteps = L,
# step = tune_param_final)
#
# par.mat <- hmc_curr$final
# accpt.par.mat.all[iter] <- hmc_curr$acc
# par_mat = par.mat
# lpd_grad_model_indep(data.rad, par_mat,
# obs_group, n.clus)
#
#
# lpd_grad_model_indep_1comp(data.rad, par_mat[, 1],
# obs_group[[1]], n.clus[1])
# lpd_grad_model_indep_1comp(data.rad, par_mat[, 2],
# obs_group[[2]], n.clus[2])
#
#
# lpd_model_indep(data.rad, par_mat,
# obs_group, n.clus)
#
# lpd_model_indep_1comp(data.rad, par_mat[, 1],
# obs_group[[1]], n.clus[1])
#
# browser()
for (j in 1:ncomp) {
# browser()
hmc_curr <-
bounded_hmc_biv(
lpr_grad =
function(par_vec_lscale) {
lpd_grad_model_indep_1comp(
data.rad, par_vec_lscale,
obs_group[[j]], n.clus[j]
)
},
init = par.mat_lscale[, j, drop = FALSE],
lower = par_lower_lscale[, j, drop = FALSE],
upper = par_upper_lscale[, j, drop = FALSE],
dep_cov = dep_cov,
dep_cov_type = dep_cov_type,
zero_cov = zero_cov,
nsteps = L,
step = tune_param_final[, j, drop = FALSE]
)
par.mat[, j] <- signif_(find_expscale(hmc_curr$final))
par.mat_lscale[, j] <- signif_(hmc_curr$final)
accpt.par.mat.all[j, iter] <- hmc_curr$acc
}
} else {
# rwmh_curr <- bounded_rwmh_biv(lpr =
# function(par_mat)
# lpd_model_indep(data.rad, par_mat,
# obs_group, n.clus),
# init = par.mat,
# lower = par_lower,
# upper = par_upper,
# dep_cov = dep_cov,
# dep_cov_type = dep_cov_type,
# zero_cov = zero_cov,
# step = tune_param)
#
# par.mat <- rwmh_curr$final
# accpt.par.mat.all[iter] <- rwmh_curr$accpt
for (j in 1:ncomp) {
rwmh_curr <- bounded_rwmh_biv(
lpr =
function(par_vec_lscale) {
lpd_model_indep_1comp(
data.rad, par_vec_lscale,
obs_group[[j]], n.clus[j]
)
},
init = par.mat_lscale[, j, drop = FALSE],
lower = par_lower_lscale[, j, drop = FALSE],
upper = par_upper_lscale[, j, drop = FALSE],
dep_cov = dep_cov,
dep_cov_type = dep_cov_type,
zero_cov = zero_cov,
step = tune_param[, j]
)
par.mat[, j] <- signif_(find_expscale(rwmh_curr$final))
par.mat_lscale[, j] <- signif_(rwmh_curr$final)
accpt.par.mat.all[j, iter] <- rwmh_curr$accpt
}
}
lprior.all[iter] <-
-0.5 * sum(par.mat_lscale[1:3, ]^2 / norm.var) +
sum(pmix.alpha * log(pi.mix))
}
# if type == "uni"
else {
if (method == "hmc") {
# hmc_curr <-
# bounded_hmc_uni(lpr_grad =
# function(par_mat)
# lpd_grad_model_indep(data.rad, par_mat,
# obs_group, n.clus),
# init = par.mat,
# lower = par_lower,
# upper = par_upper,
# nsteps = L,
# step = tune_param_final)
#
# par.mat <- hmc_curr$final
# accpt.par.mat.all[iter] <- hmc_curr$acc
for (j in 1:ncomp) {
hmc_curr <-
bounded_hmc_uni(
lpr_grad =
function(par_vec_lscale) {
lpd_grad_model_indep_1comp(
data.rad, par_vec_lscale,
obs_group[[j]], n.clus[j]
)
},
init = par.mat_lscale[, j, drop = FALSE],
lower = par_lower_lscale[, j, drop = FALSE],
upper = par_upper_lscale[, j, drop = FALSE],
nsteps = L,
step = tune_param_final[, j, drop = FALSE]
)
par.mat[, j] <- signif_(find_expscale_uni(hmc_curr$final))
par.mat_lscale[, j] <- signif_(hmc_curr$final)
accpt.par.mat.all[j, iter] <- hmc_curr$acc
}
} else {
# rwmh_curr <- bounded_rwmh_uni(lpr =
# function(par_mat)
# lpd_model_indep(data.rad, par_mat,
# obs_group, n.clus),
# init = par.mat,
# lower = par_lower,
# upper = par_upper,
# step = tune_param)
#
# par.mat <- rwmh_curr$final
# accpt.par.mat.all[iter] <- rwmh_curr$accpt
for (j in 1:ncomp) {
rwmh_curr <- bounded_rwmh_uni(
lpr =
function(par_vec_lscale) {
lpd_model_indep_1comp(
data.rad, par_vec_lscale,
obs_group[[j]], n.clus[j]
)
},
init = par.mat_lscale[, j, drop = FALSE],
lower = par_lower_lscale[, j, drop = FALSE],
upper = par_upper_lscale[, j, drop = FALSE],
step = tune_param[, j]
)
par.mat[, j] <- signif_(find_expscale_uni(rwmh_curr$final))
par.mat_lscale[, j] <- signif_(rwmh_curr$final)
accpt.par.mat.all[j, iter] <- rwmh_curr$accpt
}
}
lprior.all[iter] <-
-0.5 * sum(par.mat_lscale[1, ]^2 / norm.var) +
sum(pmix.alpha * log(pi.mix))
}
# do permutation sampling only after tuning
if (perm_sampling & iter > n.burnin) {
rand_perm <- sample(1:ncomp)
clus.ind <- rand_perm[clus.ind] # random label switch
post.wt <- post.wt[, rand_perm, drop = FALSE]
par.mat <- par.mat[, rand_perm, drop = FALSE]
par.mat_lscale <- par.mat_lscale[, rand_perm, drop = FALSE]
pi.mix <- pi.mix[rand_perm, drop = FALSE]
par.mat.all.order[, , iter] <- par.mat[, order(rand_perm), drop = FALSE]
# needed for tuning
# if (method == "hmc")
tune_param <- tune_param[, rand_perm, drop = FALSE]
}
# browser()
llik_contri <- llik_model_contri(data.rad, par.mat, pi.mix)
if (return_llik_contri) {
llik_contri_all[, iter] <- llik_contri
}
llik.all[iter] <- sum(llik_contri)
lpd.all[iter] <- llik.all[iter] + lprior.all[iter]
par.mat.all[, , iter] <- par.mat
par.mat_lscale.all[, , iter] <- par.mat_lscale
pi.mix.all[, iter] <- pi.mix
mem_prob_all[, , iter] <- post.wt
clus_ind.all[, iter] <- clus.ind
tune_param_all[, iter] <- c(tune_param)
# tuning tune_param during burnin
if (autotune & iter >= nterms & iter %% 10 == 0 &
iter <= iter.tune) {
ave_accpt_all <- rep(0, ncomp)
for (j in 1:ncomp) {
ave_accpt_all[j] <-
ave_accpt <- sum(accpt.par.mat.all[j, (iter - nterms + 1):iter]) / nterms
if (ave_accpt > accpt.prob.upper) {
tune_param[, j] <- tune_param[, j] * (1 + tune.incr)
ntunes_up[j] <- ntunes_up[j] + 1
tune_status[j, iter] <- 1
} else if (ave_accpt < accpt.prob.lower) {
tune_param[, j] <- tune_param[, j] * (1 - tune.incr)
ntunes_down[j] <- ntunes_down[j] + 1
tune_status[j, iter] <- -1
}
if (iter %% nterms == 0) {
par.sd <- apply(par.mat_lscale.all[, j, (iter - nterms + 1):iter], 1, sd)
mean_par.sd <- sum(par.sd) / (npar_1_comp)
mean_tune_param_j <- sum(tune_param[, j]) / (npar_1_comp)
if (mean_par.sd > 0) {
tune_param[, j] <- 0.5 * par.sd / mean_par.sd * mean_tune_param_j + 0.5 * tune_param[, j]
}
# else {
# par.sd <- apply(par.mat.all.order[, , (iter-nterms+1):iter], c(1:2), sd)
# par.sd <- par.sd[, rand_perm]
# mean_par.sd <- sum(par.sd)/(npar_1_comp*ncomp)
# mean_tune_param <- sum(tune_param)/(npar_1_comp*ncomp)
# if(mean_par.sd > 0)
# tune_param <- 0.5*par.sd/mean_par.sd*mean_tune_param + 0.5*tune_param
#
# cat(ave_accpt, tune_param[1,],lpd.all[iter], "\n")
# }
}
}
# if(iter %% nterms == 0)
# cat(paste("(", paste0(ave_accpt_all, collapse = ","), ")"),
# tune_param[1,], lpd.all[iter], "\n")
}
# if (method == "hmc" & autotune & iter >= nterms
# & iter %% 5 == 0 & iter <= n.burnin) {
# acr <- cor(lpd.all[(iter-nterms+2):iter],
# lpd.all[(iter-nterms+1):(iter-1)])
# if (acr < acr.lower) {
# L <- ceiling(L*(1 - tune.incr))
# } else if (acr > acr.upper) {
# L <- ceiling(L*(1 + tune.incr))
# }
# }
if (show.progress) {
message <- paste0(
"Progress: ", round(iter / n.iter * 100), "% ",
ifelse(iter <= n.burnin,
"(Burn-in)",
"(Sampling)"
)
)
if (progress.backend == "tcltk") {
tcltk::setTkProgressBar(pb, iter, label = message)
} else if (progress.backend == "txt") {
utils::setTxtProgressBar(pb, value = iter, label = message)
}
}
}
if (grepl(method, "hmc")) {
epsilon_ave <- epsilon <- mean(tune_param_all)
if (L.random) {
L_ave <- sum(L_vec) / n.iter
} else {
L_ave <- L
}
}
if (grepl(method, "rwmh")) {
propscale_final <- propscale <- rowSums(tune_param_all) / n.iter
}
allpar_val <- array(1, dim = c(npar_1_comp + 1, ncomp, n.iter))
allpar_val[1, , ] <- pi.mix.all
allpar_val[-1, , ] <- par.mat.all
rm(pi.mix.all, par.mat.all)
allpar_name <- c("pmix", modelpar.names)
dimnames(allpar_val)[[1]] <- c("pmix", modelpar.names)
if (!return_tune_param) {
rm(tune_param_all)
tune_param_all <- NULL
}
result <- list(
"par.value" = allpar_val, # [, , final_iter_set],
"par.name" = allpar_name,
"llik.contri" = llik_contri_all, # [, final_iter_set],
"llik" = llik.all, # [final_iter_set],
"lpd" = lpd.all, # [final_iter_set],
"lprior" = lprior.all, # [final_iter_set],
"accpt.modelpar" = accpt.par.mat.all, # [final_iter_set],
"clus.ind" = clus_ind.all, # [, final_iter_set],
"mem.prob" = mem_prob_all, # [, , final_iter_set],
"epsilon" = epsilon_ave,
"L" = L_ave,
"propscale" = propscale_final,
"tune_param" = tune_param_all,
"par.upper" = par_upper,
"par.lower" = par_lower,
"tcltk_fail" = tcltk_fail
)
# if (show.progress & progress.backend == "tcltk")
if (show.progress) close(pb)
result
}
# seed <- floor(runif(1, 1, 100))
# browser()
# generate three chains in parallel, if possible
if (chains_parallel) {
res_list <- future_lapply(1:n.chains,
function(ii) {
run_MC(
starting[[ii]],
L[ii], ii
)
},
future.seed = TRUE
)
} else {
res_list <- lapply(
1:n.chains,
function(ii) run_MC(starting[[ii]], L[ii], ii)
)
}
# if (res_list[[1]]$tcltk_fail) {
# msg <- paste(
# "tcltk could not be loaded; 'show.progress' was set to FALSE.",
# "Consider setting progress.backend = \"txt\" for a text-based",
# "progress bar. See ?fit_angmix for more details."
# )
# warning()
# }
# combine the results from the lists
allpar_val <- array(0, dim = c(npar_1_comp + 1, ncomp, n.iter, n.chains))
llik_all <- lprior_all <- lpd_all <- matrix(0, n.iter, n.chains)
accpt.modelpar_all <- array(0, dim = c(ncomp, n.iter, n.chains))
clus.ind_all <- array(0, dim = c(n.data, n.iter, n.chains))
if (return_llik_contri) {
llik_contri_all <- array(0, dim = c(n.data, n.iter, n.chains))
} else {
llik_contri_all <- NULL
}
mem_prob_all <- array(1, dim = c(n.data, ncomp, n.iter, n.chains))
if (return_tune_param) {
tune_param_all <- array(NA, dim = c(length(tune_param), n.iter, n.chains))
} else {
tune_param_all <- NULL
}
for (j in 1:n.chains) {
allpar_val[, , , j] <- res_list[[j]]$par.value
llik_all[, j] <- res_list[[j]]$llik
lprior_all[, j] <- res_list[[j]]$lprior
lpd_all[, j] <- res_list[[j]]$lpd
accpt.modelpar_all[, , j] <- res_list[[j]]$accpt.modelpar
clus.ind_all[, , j] <- res_list[[j]]$clus.ind
if (return_llik_contri) {
llik_contri_all[, , j] <- res_list[[j]]$llik.contri
}
mem_prob_all[, , , j] <- res_list[[j]]$mem.prob
if (return_tune_param) {
tune_param_all[, , j] <- res_list[[j]]$tune_param
}
}
if (method == "hmc") {
epsilon_final <- do.call(cbind, lapply(res_list, function(x) x$epsilon))
L_final <- unlist(lapply(res_list, function(x) x$L))
propscale_final <- NULL
} else {
propscale_final <- do.call(cbind, lapply(res_list, function(x) x$propscale))
epsilon_final <- NULL
L_final <- NULL
}
out <- list(
"par.value" = allpar_val,
"clus.ind" = clus.ind_all,
"par.name" = res_list[[1]]$par.name,
"modelpar.lower" = res_list[[1]]$par.lower,
"modelpar.upper" = res_list[[1]]$par.upper,
"return_llik_contri" = return_llik_contri,
"llik.contri" = llik_contri_all,
"mem.prob" = mem_prob_all,
"llik" = llik_all,
"lpd" = lpd_all,
"lprior" = lprior_all,
"accpt.modelpar" = accpt.modelpar_all,
"model" = curr.model,
"method" = method,
"perm_sampling" = perm_sampling,
"epsilon.random" = epsilon.random,
"epsilon" = epsilon_final,
"L.random" = L.random,
"L" = L_final,
"iter.tune" = iter.tune,
"propscale" = propscale_final,
"tune_param" = tune_param_all,
"return_tune_param" = return_tune_param,
"type" = type,
"data" = data.rad,
"cov.restrict" = cov.restrict,
"pmix.alpha" = pmix.alpha,
"norm.var" = norm.var,
"n.data" = n.data,
"ncomp" = ncomp,
"n.chains" = n.chains,
"n.iter" = n.iter,
"n.burnin" = n.burnin,
"thin" = thin,
"n.iter.final" = n.iter.final,
"final_iter" = final_iter_set,
"int.displ" = int.displ,
"qrnd_grid" = qrnd_grid,
"omega.2pi" = omega.2pi
)
class(out) <- "angmcmc"
out
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.