Nothing
#' MH sampler for BDP-DW method
#'
#' Obtain samples of the posterior of the Whittle likelihood in conjunction with a Bernstein-Dirichlet prior on the spectral density.
#' @details Further detail can be found in the simulation study section in the references papers.
#' @param data numeric vector.
#' @param m window size needed to calculate moving periodogram.
#' @param likelihood_thinning the thinning factor of the dynamic Whittle likelihood.
#' @param mcmc_params a list generated by \link{bdp_dw_mcmc_params_gen}.
#' @param prior_params a list generated by \link{bdp_dw_prior_params_gen}.
#' @param monitor a Boolean value (default FALSE) indicating whether to display the real-time status
#' @param print_interval If monitor = TRUE, then this value indicates the number of iterations after which a status is printed to console;
#' If monitor = FALSE, it does not have any effect
#' @return list containing the following fields:
#'
#' \item{k1,k2,tau,V,W1,W2}{posterior traces of PSD parameters}
#' \item{tim}{total run time}
#' \item{prior_params}{the specifications of the prior}
#' @references Y. Tang et al. (2023)
#' \emph{Bayesian nonparametric spectral analysis of locally stationary processes}
#' ArXiv preprint
#' <arXiv:2303.11561>
#'
#' @importFrom Rcpp evalCpp sourceCpp
#' @importFrom stats cov pbeta qlogis rbinom rnorm rpois runif
#' @useDynLib beyondWhittle, .registration = TRUE
#' @keywords internal
bdp_dw_mcmc <- function(data,
m,
likelihood_thinning = 1,
mcmc_params,
prior_params,
monitor = FALSE,
print_interval = 100)
{
# MCMC parameters
stopifnot(!is.null(mcmc_params$Ntotal)); stopifnot(mcmc_params$Ntotal>0)
Ntotal <- mcmc_params$Ntotal
stopifnot(!is.null(mcmc_params$burnin)); stopifnot(mcmc_params$burnin>=0 && mcmc_params$burnin<Ntotal)
burnin <- mcmc_params$burnin
stopifnot(!is.null(mcmc_params$thin)); stopifnot(mcmc_params$thin>=1)
thin <- mcmc_params$thin
#
# Adaptive MCMC parameters
stopifnot(!is.null(mcmc_params$adaption.batchSize)); stopifnot(mcmc_params$adaption.batchSize>0)
adaption.batchSize <- mcmc_params$adaption.batchSize
stopifnot(!is.null(mcmc_params$adaption.targetAcceptanceRate)); stopifnot(mcmc_params$adaption.targetAcceptanceRate>0 && mcmc_params$adaption.targetAcceptanceRate<1)
adaption.targetAcceptanceRate <- mcmc_params$adaption.targetAcceptanceRate
# PRIOR PAREMETERS
stopifnot(!is.null(prior_params$M)); stopifnot(prior_params$M > 0)
M <- prior_params$M
stopifnot(!is.null(prior_params$g0.alpha) && !is.null(prior_params$g0.beta)); stopifnot(prior_params$g0.alpha > 0 && prior_params$g0.beta > 0)
g0.alpha <- prior_params$g0.alpha
g0.beta <- prior_params$g0.beta
stopifnot(!is.null(prior_params$k1.theta)); stopifnot(prior_params$k1.theta > 0)
k1.theta <- prior_params$k1.theta
stopifnot(!is.null(prior_params$k2.theta)); stopifnot(prior_params$k2.theta > 0)
k2.theta <- prior_params$k2.theta
stopifnot(!is.null(prior_params$tau.alpha) && !is.null(prior_params$tau.beta)); stopifnot(prior_params$tau.alpha > 0 && prior_params$tau.beta > 0)
tau.alpha <- prior_params$tau.alpha
tau.beta <- prior_params$tau.beta
stopifnot(!is.null(prior_params$k1max)); stopifnot(prior_params$k1max > 0)
k1max <- prior_params$k1max
stopifnot(!is.null(prior_params$k2max)); stopifnot(prior_params$k2max > 0)
k2max <- prior_params$k2max
stopifnot(!is.null(prior_params$L)); stopifnot(prior_params$L > 0)
L <- prior_params$L
stopifnot(!is.null(prior_params$bernstein1_l) && !is.null(prior_params$bernstein1_r)); stopifnot(prior_params$bernstein1_l >= 0 && prior_params$bernstein1_r <= 1)
bernstein1_l <- prior_params$bernstein1_l
bernstein1_r <- prior_params$bernstein1_r
stopifnot(!is.null(prior_params$bernstein2_l) && !is.null(prior_params$bernstein2_r)); stopifnot(prior_params$bernstein2_l >= 0 && prior_params$bernstein2_r <= 1)
bernstein2_l <- prior_params$bernstein2_l
bernstein2_r <- prior_params$bernstein2_r
#
stopifnot(is.logical(monitor))
stopifnot(!is.null(print_interval)); stopifnot(print_interval>0)
stopifnot(!is.null(likelihood_thinning)); stopifnot(likelihood_thinning %in% 1:3)
#
n <- length(data) - 2*m
# moving Fourier coefs and time_grid
mf_list <- local_moving_FT_zigzag(data, m, likelihood_thinning)
time_grid <- mf_list$time_grid # original scale
# moving periodograms (with scale 1/(2*pi))
FZ <- abs(mf_list$MF)^2/(2*pi)
# grid points
u1 <- time_grid/n # re-scaled
u2 <- (2 * (1:m))/(2*m+1)
# Bernstein basis
db.list_1 <- dbList_dw_Bern(u1, k1max, bernstein1_l, bernstein1_r)
db.list_2 <- dbList_dw_Bern_for_lambda(u2, k2max, bernstein2_l, bernstein2_r, m, time_grid)
# Open objects for storage
tau <- rep(NA, Ntotal)
tilde.E <- matrix(NA, nrow = L + 1, ncol = Ntotal) # tilde.E = log(E)
W1 <- matrix(NA, nrow = L, ncol = Ntotal)
W2 <- matrix(NA, nrow = L, ncol = Ntotal)
k1 <- rep(NA, Ntotal)
k2 <- rep(NA, Ntotal)
lpost_storage <- rep(NA, Ntotal)
# Starting values
tau[1] <- tau.hat <- mean(FZ)
k1[1] <- round(k1max / 2)
k2[1] <- round(k2max / 2)
beta_basis_1_k <- db.list_1[[k1[1]]]
beta_basis_2_k <- db.list_2[[k2[1]]]
tilde.E[, 1] <- log(rep(1/(L+1), L+1)) # log transformation s.t. the domain becomes the whole space
W1[, 1] <- seq(from=1/(2*k1[1]), to=1-1/(2*k1[1]), length.out=L)
W2[, 1] <- seq(from=1/(2*k2[1]), to=1-1/(2*k2[1]), length.out=L)
p.current <- qgamma(1 - cumsum(exp(tilde.E[-(L+1), 1]))/sum(exp(tilde.E[, 1])), shape = M/L, rate = 1)
selector1.current <- findInterval(W1[, 1], (1:(k1[1]))/(k1[1]), left.open = T) + 1 # amount to ceiling(k1 * w1) but safer
# mat1.current <- beta_basis_1_k[selector1.current, ]
selector2.current <- findInterval(W2[, 1], (1:(k2[1]))/(k2[1]), left.open = T) + 1 # amount to ceiling(k2 * w2) but safer
# mat2.current <- beta_basis_2_k[selector2.current, ]
lp.current <- lprior_dw(tilde.E[,1], W1[, 1], W2[, 1], k1[1], k2[1], tau[1],
g0.alpha, g0.beta, k1.theta, k2.theta, tau.alpha, tau.beta)
npsd.current <- qpsd_cal_cpp_expedited(beta_basis_1_k,
beta_basis_2_k,
p.current,
integer(ncol(beta_basis_1_k)),
selector1.current - 1,
selector2.current - 1)/sum(p.current)
ll.current <- llike_dw(FZ, npsd.current, tau[1])
lpost_storage[1] <- lp.current + ll.current
# proposal variances
# initial Metropolis proposal parameters for W1, W2
eps_W1 <- eps_W2 <- seq(1, L) / (seq(1, L) + 2 * sqrt(n))
eps_E <- seq(1, L+1) / (seq(1, L+1) + 2 * sqrt(n))
eps_tau <- 5
lsd_E <- log(eps_E) / 2
lsd_W1 <- log(eps_W1) / 2
lsd_W2 <- log(eps_W2) / 2
####
# Metropolis-within-Gibbs sampler
tim.start <- Sys.time()
tim0 <- proc.time()
for (i in 1:(Ntotal-1)) {
if (monitor){
if ((i+1)%%print_interval == 0) {
tim <- proc.time()
print_mcmc_state(i+1, Ntotal, tim, tim0)
}
}
# adaptation
if ((i > 1) && (i <= burnin) && (i %% adaption.batchSize == 1)) {
batch <- (i - 1):(i - adaption.batchSize)
adaption.delta <- min(0.05, 1/sqrt(i)) # c.f. Rosenthal
### tilde.E
batch.E <- tilde.E[, batch]
batch.E.acceptanceRate <- apply(batch.E, 1, acceptanceRate)
lsd_E <- lsd_E + ((batch.E.acceptanceRate > adaption.targetAcceptanceRate)*2-1) * adaption.delta
eps_E <- exp(2*lsd_E)
### W1
batch.W1 <- W1[, batch]
batch.W1.acceptanceRate <- apply(batch.W1, 1, acceptanceRate)
lsd_W1 <- lsd_W1 + ((batch.W1.acceptanceRate > adaption.targetAcceptanceRate)*2-1) * adaption.delta
eps_W1 <- exp(2*lsd_W1)
### W2
batch.W2 <- W2[, batch]
batch.W2.acceptanceRate <- apply(batch.W2, 1, acceptanceRate)
lsd_W2 <- lsd_W2 + ((batch.W2.acceptanceRate > adaption.targetAcceptanceRate)*2-1) * adaption.delta
eps_W2 <- exp(2*lsd_W2)
}
# STEP 1.1: Sample Bernstein polynomial degree k1
k1.star <- k1[i] + (rpois(1, lambda = 1) + 1) * (2*rbinom(1, 1, 0.5) - 1)
k1.star <- (k1.star-1)%%as.integer(k1max) + 1 # circular symmetric proposal
beta_basis_1_k_star <- db.list_1[[k1.star]]
selector1.star <- findInterval(W1[, i], (1:k1.star)/k1.star, left.open = T) + 1
lp.k1.star <- lprior_dw(tilde.E[,i], W1[, i], W2[, i], k1.star, k2[i], tau[i],
g0.alpha, g0.beta, k1.theta, k2.theta, tau.alpha, tau.beta)
npsd.k1.star <- qpsd_cal_cpp_expedited(beta_basis_1_k_star,
beta_basis_2_k,
p.current,
integer(ncol(beta_basis_1_k_star)),
selector1.star - 1,
selector2.current - 1)/sum(p.current)
ll.k1.star <- llike_dw(FZ, npsd.k1.star, tau[i])
#####
# Accept/reject
alpha_k1 <- min(0, (lp.k1.star + ll.k1.star) - (lp.current + ll.current))
if (isTRUE(log(runif(1, 0, 1)) < alpha_k1)) {
k1[i + 1] <- k1.star # Accept k1.star
lp.current <- lp.k1.star
npsd.current <- npsd.k1.star
ll.current <- ll.k1.star
beta_basis_1_k <- beta_basis_1_k_star
selector1.current <- selector1.star
} else {
k1[i + 1] <- k1[i] # Reject and use previous
}
# STEP 1.2: Sample Bernstein polynomial degree k2
k2.star <- k2[i] + (rpois(1, lambda = 1) + 1) * (2*rbinom(1, 1, 0.5) - 1)
k2.star <- (k2.star-1)%%as.integer(k2max) + 1 # circular symmetric proposal
beta_basis_2_k_star <- db.list_2[[k2.star]]
selector2.star <- findInterval(W2[, i], (1:k2.star)/k2.star, left.open = T) + 1
lp.k2.star <- lprior_dw(tilde.E[,i], W1[, i], W2[, i], k1[i+1], k2.star, tau[i],
g0.alpha, g0.beta, k1.theta, k2.theta, tau.alpha, tau.beta)
npsd.k2.star <- qpsd_cal_cpp_expedited(beta_basis_1_k,
beta_basis_2_k_star,
p.current,
integer(ncol(beta_basis_1_k)),
selector1.current - 1,
selector2.star - 1)/sum(p.current)
ll.k2.star <- llike_dw(FZ, npsd.k2.star, tau[i])
#####
# Accept/reject
alpha_k2 <- min(0, (lp.k2.star + ll.k2.star) - (lp.current + ll.current))
if (isTRUE(log(runif(1, 0, 1)) < alpha_k2)) {
k2[i + 1] <- k2.star # Accept k.star
lp.current <- lp.k2.star
npsd.current <- npsd.k2.star
ll.current <- ll.k2.star
beta_basis_2_k <- beta_basis_2_k_star
selector2.current <- selector2.star
} else {
k2[i + 1] <- k2[i] # Reject and use previous
}
# can be implemented using cpp armadillo
mat12 <- beta_basis_1_k[selector1.current, ] * beta_basis_2_k[selector2.current, ]
# sampling tilde.E
tilde.E.star <- tilde.E[, i]
for(l in 1:(L+1)){
if (l > 1) {
tilde.E.star[l-1] <- tilde.E[l-1, i + 1]
tilde.E.star[l-1] <- tilde.E[l-1, i + 1]
}
tilde.E.star[l] <- tilde.E.star[l] + eps_E[l] * rnorm(1)
lp.E.star <- lprior_dw(tilde.E.star, W1[,i], W2[,i], k1[i+1], k2[i+1], tau[i],
g0.alpha, g0.beta, k1.theta, k2.theta, tau.alpha, tau.beta)
p.star <- qgamma(1 - cumsum(exp(tilde.E.star[-(L+1)]))/sum(exp(tilde.E.star)), shape = M/L, rate = 1)
npsd.E.star <- colSums(p.star * mat12)/sum(p.star)
ll.E.star <- llike_dw(FZ, npsd.E.star, tau[i])
#####
# Accept/reject
alpha_E <- min(0, (ll.E.star - ll.current) + (lp.E.star - lp.current))
if (isTRUE(log(runif(1, 0, 1)) < alpha_E)) {# accept the proposed value
tilde.E[l, i+1] <- tilde.E.star[l]
lp.current <- lp.E.star
p.current <- p.star
ll.current <- ll.E.star
npsd.current <- npsd.E.star
} else {# Reject and use previous
tilde.E[l, i+1] <- tilde.E[l, i]
}
}
unip.current <- p.current/sum(p.current)
unipmat12 <- unip.current * mat12
# npsd.current <- colSums(unipmat12)
unipmat12.star <- unipmat12
# updating W1 and W2
W1.star <- W1[, i]
W2.star <- W2[, i]
for (l in 1:L) {
if (l > 1) {
W1.star[l-1] <- W1[l-1, i + 1]
W2.star[l-1] <- W2[l-1, i + 1]
}
###---first W1, then W2---------------
# Uniform proposal (W1[,i] - eps, W1[,i] + eps) on (0,1)
W1.star[l] <- W1.star[l] + eps_W1[l] * rnorm(1)
W1.star[l] <- W1.star[l] - floor(W1.star[l])
# W1.star[l] <- runif(1, W1.star[l] - eps_W1[l], W1.star[l] + eps_W1[l])
# W1.star[l] <- W1.star[l] - floor(W1.star[l])
selector1l.star <- findInterval(W1.star[l], (1:(k1[i+1]))/(k1[i+1]), left.open = T) + 1
mat12lvec.star <- beta_basis_1_k[selector1l.star, ] * beta_basis_2_k[selector2.current[l], ]
# npsd.W.star <- npsd.current - unipmat12[l, ] + unip.current[l] * mat12lvec.star
unipmat12.star[l, ] <- unip.current[l] * mat12lvec.star
# npsd.W.star <- colSums(unipmat12.star)
npsd.W.star <- npsd.current - unipmat12[l, ] + unipmat12.star[l, ]
lp.W.star <- lprior_dw(tilde.E[,i+1], W1.star, W2.star, k1[i+1], k2[i+1], tau[i],
g0.alpha, g0.beta, k1.theta, k2.theta, tau.alpha, tau.beta)
ll.W.star <- llike_dw(FZ, npsd.W.star, tau[i])
#####
# Accept/reject
alpha_W <- min(0, (ll.W.star - ll.current) + (lp.W.star - lp.current))
if (isTRUE(log(runif(1, 0, 1)) < alpha_W)) {# accept the proposed value
W1[l, i+1] <- W1.star[l]
lp.current <- lp.W.star
ll.current <- ll.W.star
npsd.current <- npsd.W.star
unipmat12[l, ] <- unipmat12.star[l, ]
selector1.current[l] <- selector1l.star
} else {# Reject and use previous
W1[l, i+1] <- W1[l, i]
unipmat12.star[l, ] <- unipmat12[l, ]
}
##--------------------------
# Uniform proposal (W1[,i] - eps, W1[,i] + eps) on (0,1)
W2.star[l] <- W2.star[l] + eps_W2[l] * rnorm(1)
W2.star[l] <- W2.star[l] - floor(W2.star[l])
# W2.star[l] <- runif(1, W2.star[l] - eps_W2[l], W2.star[l] + eps_W2[l])
# W2.star[l] <- W2.star[l] - floor(W2.star[l])
selector2l.star <- findInterval(W2.star[l], (1:(k2[i+1]))/(k2[i+1]), left.open = T) + 1
mat12lvec.star <- beta_basis_1_k[selector1.current[l], ] * beta_basis_2_k[selector2l.star, ]
# npsd.W.star <- npsd.current - unipmat12[l, ] + unip.current[l] * mat12lvec.star
unipmat12.star[l, ] <- unip.current[l] * mat12lvec.star
# npsd.W.star <- colSums(unipmat12.star)
npsd.W.star <- npsd.current - unipmat12[l, ] + unipmat12.star[l, ]
lp.W.star <- lprior_dw(tilde.E[,i+1], W1.star, W2.star, k1[i+1], k2[i+1], tau[i],
g0.alpha, g0.beta, k1.theta, k2.theta, tau.alpha, tau.beta)
ll.W.star <- llike_dw(FZ, npsd.W.star, tau[i])
#####
# Accept/reject
alpha_W <- min(0, (ll.W.star - ll.current) + (lp.W.star - lp.current))
if (isTRUE(log(runif(1, 0, 1)) < alpha_W)) {# accept the proposed value
W2[l, i+1] <- W2.star[l]
lp.current <- lp.W.star
ll.current <- ll.W.star
npsd.current <- npsd.W.star
unipmat12[l, ] <- unipmat12.star[l, ]
selector2.current[l] <- selector2l.star
} else {# Reject and use previous
W2[l, i+1] <- W2[l, i]
unipmat12.star[l, ] <- unipmat12[l, ]
}
}
# sampling tau
tau.star <- tau.hat * runif(1, min = 0, max = eps_tau)
lp.tau.star <- lprior_dw(tilde.E[, i+1], W1[, i+1], W2[, i+1], k1[i+1], k2[i+1], tau.star,
g0.alpha, g0.beta, k1.theta, k2.theta, tau.alpha, tau.beta)
ll.tau.star <- llike_dw(FZ, npsd.current, tau.star)
#####
# Accept/reject
alpha_tau <- min(0, (ll.tau.star - ll.current) + (lp.tau.star - lp.current))
if (isTRUE(log(runif(1, 0, 1)) < alpha_tau)) {# Accept tau.star
tau[i + 1] <- tau.star
lp.current <- lp.tau.star
ll.current <- ll.tau.star
lpost_storage[i+1] <- lp.current + ll.current # tracing log posterior
} else {# Reject and use previous
tau[i + 1] <- tau[i]
lpost_storage[i+1] <- lpost_storage[i] # tracing log posterior
}
} # END: MCMC loop
tim.end <- Sys.time()
# Which iterations to keep
keep <- seq(burnin + 1, Ntotal, by = thin)
k1 <- k1[keep]
k2 <- k2[keep]
tau <- tau[keep]
E <- exp(tilde.E[, keep])
W1 <- W1[, keep]
W2 <- W2[, keep]
lpost_storage <- lpost_storage[keep]
return(list(k1 = k1, k2 = k2, tau = tau, E = E, W1 = W1, W2 = W2,
tim = tim.end - tim.start, prior_params = prior_params, lpost = lpost_storage, data = data)
)
}
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.