R/dynamicWhittle_MCMC_zigzag.R

Defines functions bdp_dw_mcmc

Documented in bdp_dw_mcmc

#' 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)
         )
}

Try the beyondWhittle package in your browser

Any scripts or data that you put into this service are public.

beyondWhittle documentation built on April 3, 2025, 7:49 p.m.