runmcmc_cp1_left {bulletcp}

Estimate a posterior distribution of data conditional on a left groove and no right groove.

Description

This function runs a random walk metropolis within Gibbs algorithm to estimate the posterior distribution of the value of the changepoint as well as the parameters fit in each multivariate normal distribution on either side of the changepoint. The covariance matrices are both based on the exponential covariance function. This functions assumes equally spaced locations ("x" values in the "data" argument). The distribution to the left of the changepoint has a mean that is a linear function of the distance from the center of the data.

Usage

runmcmc_cp1_left(data, iter, start.vals, prop_var, cp_prop_var,
  tol_edge = 50, warmup = 500, verbose = FALSE)

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 interations after warmup.

start.vals

List with elements "sigma", "l", "cp", "beta", and "intercept." "sigma" and "l" are 2 element vectors where the first element is for the data to the left of the changepoint. "cp" is the changepoint starting value. "beta" and "intercept" are the slope and intercept starting values for the mean of the data model to the left of the changepoint. which parameterize the covariance matrix.

prop_var

A two element list of the proposal variance-covariance matrices for the random walk metropolis algorithm(s). The first element is for the data to the left of the changepoint.

cp_prop_var

The proposal variance for the changepoint.

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 the proposal is within a distance of 10 x-values from either edge.

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.

verbose

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

Value

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. "gp_prop_var" and "cp_prop_var" are the tuned proposal variances for the metropolis steps.

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 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
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)))

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

# define the proposal variance for the RWMH step sampling the changepoint
cp_prop_var <- 10^2

# run Gibbs MCMC for the left GEA model
set.seed(1111)
m1cp_left <- runmcmc_cp1_left(data = fake_groove, iter = 500, warmup = 100,
                           start.vals = start.vals$left,
                           prop_var = prop_var$left,
                           cp_prop_var = cp_prop_var,
                           verbose = FALSE, tol_edge = 50)