bulletcp (version 1.0.0)

runmcmc_cpall: Estimate posterior distributions for the 0, 1, or 2 changepoint case.

Description

This function runs the changepoint functions designed for the cases when there are 0, 1, or 2 changepoints. It then returns a subset of the results that are returned for each function individually. This subset of results is enough to decide the likely number of shoulders, the locations of the shoulders (if they exist), as well as the posterior samples for the changepoints for minimal diagnostic use.

Usage

runmcmc_cpall(data, iter = 8000, start.vals = NA, prop_var = NA,
  cp_prop_var = NA, tol_edge = 50, tol_cp = 1000, warmup = 500,
  verbose = FALSE, prior_numcp = rep(1/4, times = 4))

Arguments

data

Data frame with columns "x" and "y." "x" is a column of the locations of the observed residual values, y.

iter

Number of iterations after warmup.

start.vals

Starting values for the changepoint algorithm. Either NA valued or a named list of lists. If list, the names of the lists should be "cp2","cp1", and "cp0". Each list posessing one of those aforementioned names is a list of starting values identical to what would be given if the changepoint algorithm were to be run with the corresponding number of specified changepoints. List with elements "sigma", "l", "cp", "beta", and "intercept." "sigma" and "l" are 3 element vectors where the first element is for the data on the left groove. The second element is for the land engraved area, and the third element is for the right groove. "cp" is the vector of changepoint starting values. "beta" and "intercept" are two element vectors of the slope and intercept for the left and right groove engraved area respectively. If NA, default starting values will be used. Note that the changepoint starting values should always be near the edges of the data.

prop_var

Either NA valued or a list of named lists. If list, the names of the lists should be "cp2","cp1", and "cp0". Each list posessing one of those aforementioned names is a list of proposal covariance matrices identical to what would be given if the changepoint algorithm were to be run with the corresponding number of specified changepoints.

cp_prop_var

The proposal variance-covariance matrix for the changepoints. Can either be NA or a named list. If list, the names of the list items should be "cp2", "cp1" where each is the appropriate proposal variance/covariance matrix for the number of changepoints.

tol_edge

This parameter controls how close changepoint proposals can be to the edge of the data before getting automatically rejected. For example, a value of 10 means that the changepoint will be automatically rejected if either of the proposal changepoints is within a distance of 10 x-values from either edge.

tol_cp

This parameter controls how close changepoint proposals can be to each other before getting automatically rejected. For example, a value of 10 means that the changepoint will be automatically rejected if either of the proposal changepoints is within a distance of 10 x-values from either each other.

warmup

The number of warmup iterations. This should be set to a very small number of iterations, as using too many iterations as warmup risks moving past the changepoints and getting stuck in a local mode. Default is set to 500.

verbose

Logical value indicating whether to print the iteration number and the parameter proposals.

prior_numcp

This is a vector with four elements giving the prior probabilities for the zero changepoint model, the one changepoint on the left model, the one changepoint on the right model, and the two changepoint model, in that order. Note that, practically, because the likelihood values are so large, only very strong priors will influence the results.

Value

A named list containing the sampled changepoint locations for both the one and two changepoint scenarios, the posterior changepoint means, the average log pdf values from the data model under each model, the maximum log probability values under each model log likelihood values, and estimates of the maximum a posteriori changepoint value under each model.

Examples

# NOT RUN {
# 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 for the changepoints
cp_start_left <- min(fake_groove$x) + 60
cp_start_right <- max(fake_groove$x) - 60

# define list of starting values for both the left and right changepoint models
cp0.start.vals <- list("sigma" = c(1), "l" = c(10))
cp1.start.vals <- list("left" = list("sigma" = c(1,1),
                              "l" = c(10,10),
                              "cp" = c(cp_start_left),
                              "beta" = c(-1),
                              "intercept" = c(0)),
                              "right" = list("sigma" = c(1,1),
                               "l" = c(10,10),
                                "cp" = c(cp_start_right),
                                "beta" = c(1),
                                "intercept" = c(0)))
cp2.start.vals <- list("sigma" = c(1,1,1),
                "l" = c(10,10,10),
                "cp" = c(cp_start_left, cp_start_right),
                "beta" = c(-2,2),
                "intercept" = c(0,0))
start.vals <- list("cp2" = cp2.start.vals, "cp1" = cp1.start.vals, "cp0" = cp0.start.vals)

# list of starting values for each of the two MH steps
# (not sampling the changepoint) for both the left and right changepoint models
prop_var_0cp <- diag(c(1/2,1/2))
prop_var_lrcp <- list("left" = list(diag(c(1/2,1/2,1/2,1/2)),
                            diag(c(1/2,1/2))),
                            "right" = list(diag(c(1/2,1/2)),
                            diag(c(1/2,1/2,1/2, 1/2))))

prop_var_2cp <- list(diag(c(1/2,1/2,1/2,1/2)),
              diag(c(1/2,1/2)),
              diag(c(1/2,1/2,1/2,1/2)))

prop_var <- list("cp2" = prop_var_2cp, "cp1" = prop_var_lrcp, "cp0" = prop_var_0cp)

# define the proposal variance for the RWMH step sampling the changepoint
cp_prop_var <- list("cp2" = diag(c(10^2, 10^2)),
                 "cp1" = 10^2)

# prior on the number of changepoints
prior_numcp <- rep(1/4, times = 4)

set.seed(1111)
cp_gibbs <- runmcmc_cpall(data = fake_groove,
                       start.vals = start.vals,
                       prop_var = prop_var,
                       cp_prop_var = cp_prop_var,
                       verbose = FALSE,
                       tol_edge = 50,
                       tol_cp = 1000,
                       iter = 300,
                       warmup = 100,
                       prior_numcp = prior_numcp)
# }