R/changepoint_gibbs_0cp.R

In bulletcp: Automatic Groove Identification via Bayesian Changepoint Detection

#'  Estimate a posterior distribution of data conditional on zero changepoints.
#'
#' This function runs a random walk Metropolis algorithm to estimate the posterior distribution
#' of a zero mean multivariate normal distribution with an covariance matrix generated by the
#' exponential covariance function. This functions assumes equally spaced locations ("x" values in the
#' "data" argument).
#' @param data Data frame with columns "x" and "y." "x" is a column of the locations of the
#' observed residual values, y.
#' @param iter Number of interations after warmup.
#' @param start.vals List with elements "sigma" and "l" for the standard deviation and length scale
#'  which parameterize the covariance matrix.
#' @param prop_var The proposal variance-covariance matrix for the random walk metropolis algorithm.
#' @param warmup The number of initial iterations which serves two purposes: the first is to allow the
#' algorithm to wander to the area of most mass, and the second is to tune the proposal variance.
#' @param verbose Logical value indicating whether to print the iteration number and the parameter proposals.
#' @return A named list. "parameters" is a list of named parameter values each of which is a vector of length
#' "iter". "accept" gives the proportion of accepted proposals after warmup. "lp" is a vector of
#' values of the log data pdf at each sampled parameter value.
#' @importFrom stats dgamma dnorm runif var
#' @examples
#' # Fake data
#' sim_groove <- function(beta = c(-0.28,0.28), a = 125)
#' {
#'     x <- seq(from = 0, to = 2158, by = 20)
#'     med <- median(x)
#'     y <- 1*(x <= a)*(beta[1]*(x - med) - beta[1]*(a - med)) +
#'     1*(x >= 2158 - a)*(beta[2]*(x - med) - beta[2]*(2158 - a - med))
#'     return(data.frame("x" = x, "y" = y))
#' }
#'
#' fake_groove <- sim_groove()
#' # define starting values
#' start.vals <- list("sigma" = c(1), "l" = c(10))
#'
#' # proposal variance for the MH step
#' prop_var <- diag(c(1/2,1/2))
#'
#' set.seed(1111)
#' m0cp <- runmcmc_cp0(data = fake_groove, iter = 500,
#'     start.vals = start.vals,
#'     prop_var = prop_var, warmup = 100, verbose = FALSE)
#' @export

runmcmc_cp0 <- function(data, iter, start.vals, prop_var, warmup = 500, verbose = FALSE)
{
  ##data is a data frame with column x and column y

  lognormal_ou_pdf <- function(x, mu, sigma, l)
  {
    n <- length(x)
    rho <- exp(-1/l)

    return(-n/2 * log(2 * pi) - n * log(sigma) - ((n - 1)/2) * log(1 - rho^2)
           - 1/2 * 1/(sigma^2 * (1 - rho^2)) * ((x[1] - mu[1])^2 + (x[n] - mu[n])^2 + (1 + rho^2) * sum((x[2:(n-1)] - mu[2:(n-1)])^2)
                                                - 2 * rho * sum((x[1:(n-1)] - mu[1:(n-1)]) * (x[2:n] - mu[2:n]))))
  }

  ## initialize parameter list
  par <- list()
  par$sigma <- matrix(nrow = warmup + 1, ncol = 1) ## the variance of the GP
  par$sigma[1,] <- start.vals$sigma

  par$l <- matrix(nrow = warmup + 1, ncol = 1) ## length scale of the GP
  par$l[1,] <- start.vals$l

  # par$tau <- matrix(nrow = warmup + 1, ncol = 1) ## nugget of the data model
  # par$tau[1,] <- start.vals$tau

  ## range on the x-axis of data
  interval <- range(data$x)

  ## current values of parameters
  sigma <- start.vals$sigma
  l <- start.vals$l
  # tau <- start.vals$tau

  ## initialize acceptance rates
  accept <- list()
  accept$gp_par <- matrix(data = c(0), nrow = 1, ncol = 1)

  ## warmup iterations
  for(i in 1:(warmup)){

  ## make proposal for MH steps for GP parameters
  prop <- as.numeric(mvtnorm::rmvnorm(n = 1, mean = c(sigma, l), sigma = prop_var))
  if(verbose == TRUE)
  {
    print(paste("iteration: ",i))
    print(paste("GP parameter proposal: ", prop))
  }

  ## if the proposal produces negative values, reject
  if(any(prop <= 0))
  {
    # par$sigma[i,j] <- sigma[j]
    # par$l[i,j] <- l[j]
    # par$tau[i,j] <- tau[j]

    ## update GP parameters
    par$sigma[i + 1,] <- sigma
    par$l[i + 1,] <- l
    # par$tau[i + 1,] <- tau
  }
  # if both proposal values are positive, compute acceptance probability
  else{
    temp_dat <- data$y

    ## acceptance probability is just the ratio of posteriors at the proposed value
    ##   to the current value
    ## proposal doesn't appear because it should cancel as it is symmetric
    log_accept_ratio <- lognormal_ou_pdf(x = temp_dat, mu = rep(0, times = length(temp_dat)), sigma = prop[1], l = prop[2]) + #log likelihood
      dgamma(x = prop[2], shape = 3, rate = 5, log = TRUE) + # log density of the length scale
      dnorm(x = prop[1], mean = 0, sd = 1, log = TRUE) - # log density of the standard deviation
      (lognormal_ou_pdf(x = temp_dat, mu = rep(0, times = length(temp_dat)), sigma = sigma, l = l) +
         dgamma(x = l, shape = 3, rate = 5, log = TRUE) +
         dnorm(x = sigma, mean = 0, sd = 1, log = TRUE))

    ## accept the proposal with the appropriate probability
    if(log(runif(n = 1, min = 0, max = 1)) <= log_accept_ratio)
    {
      sigma <- prop[1]
      l <- prop[2]
      # tau <- prop[3]
    }
    ## update GP parameters
    par$sigma[i + 1,] <- sigma
    par$l[i + 1,] <- l
    # par$tau[i + 1,] <- tau
  }
    #print(i)
  }
  ###########################################################
  ## End warmup
  ###########################################################

  ## documentation here is sparse as the procedure is the same as that during warmup

  ## tune metropolis proposal variances
  prop_var <- 2.4^2 * var(cbind(par$sigma[round(warmup/2):warmup,1], par$l[round(warmup/2):warmup,1])) / 2

  ## reinitialize parameter list
  lp <- numeric()
  lpost <- numeric() ## log posterior value
  par <- list()
  par$sigma <- matrix(nrow = iter + 1, ncol = 1) ## the variance of the GP
  par$sigma[1,] <- sigma

  par$l <- matrix(nrow = iter + 1, ncol = 1) ## length scale of the GP
  par$l[1,] <- l

  # par$tau <- matrix(nrow = iter + 1, ncol = 1) ## nugget of the data model
  # par$tau[1,] <- tau


  ## gibbs sampling iterations
  for(i in 1:(iter))
  {
    ## given changepoints make proposal for MH steps for GP parameters
    prop <- as.numeric(mvtnorm::rmvnorm(n = 1, mean = c(sigma, l), sigma = prop_var))
    if(verbose == TRUE)
    {
      print(paste("iteration: ",i))
      print(paste("GP parameter proposal: ", prop))
    }

    ## skip this chunk of data if the proposal produces negative values
    if(any(prop <= 0))
    {
      ## update GP parameters
      par$sigma[i + 1,] <- sigma
      par$l[i + 1,] <- l
      # par$tau[i + 1,] <- tau

      temp_dat <- data$y

      ## compute and save the log likelihood values to be used in model selection
      lp[i] <- lognormal_ou_pdf(x = temp_dat, mu = rep(0, times = length(temp_dat)), sigma = sigma, l = l)

    }
    else{
      temp_dat <- data$y

      ## proposal doesn't appear because it should cancel
      log_accept_ratio <- lognormal_ou_pdf(x = temp_dat, mu = rep(0, times = length(temp_dat)), sigma = prop[1], l = prop[2]) +
        dgamma(x = prop[2], shape = 3, rate = 5, log = TRUE) +
        dnorm(x = prop[1], mean = 0, sd = 1, log = TRUE) -
        (lognormal_ou_pdf(x = temp_dat, mu = rep(0, times = length(temp_dat)), sigma = sigma, l = l) +
           dgamma(x = l, shape = 3, rate = 5, log = TRUE) +
           dnorm(x = sigma, mean = 0, sd = 1, log = TRUE))

      ## compute and save the log likelihood values to be used in model selection
      lp[i] <- lognormal_ou_pdf(x = temp_dat, mu = rep(0, times = length(temp_dat)), sigma = sigma, l = l)
      lpost[i] <- lp[i] + dgamma(x = l, shape = 3, rate = 5, log = TRUE) +
        dnorm(x = sigma, mean = 0, sd = 1, log = TRUE)

      if(log(runif(n = 1, min = 0, max = 1)) <= log_accept_ratio)
      {
        accept$gp_par[1,] <- accept$gp_par[1,] + 1/iter
        sigma <- prop[1]
        l <- prop[2]
        # tau <- prop[3]

        ## compute and save the log likelihood values to be used in model selection
        lp[i] <- lognormal_ou_pdf(x = temp_dat, mu = rep(0, times = length(temp_dat)), sigma = prop[1], l = prop[2])
        lpost[i] <- lp[i] + dgamma(x = prop[2], shape = 3, rate = 5, log = TRUE) +
          dnorm(x = prop[1], mean = 0, sd = 1, log = TRUE)
      }
      ## update GP parameters
      par$sigma[i + 1,] <- sigma
      par$l[i + 1,] <- l
      # par$tau[i + 1,] <- tau
    }
  }
  return(list("parameters" = par, "accept" = accept, "lp" = lp, "lpost" = lpost))
}