# SHORT-MET.R - Functions for short-cut Metropolis and related methods. # COUNT OF NUMBER OF LOG PROBABILITY EVALUATIONS. lpr.evals <<- 0 # DO ONE METROPOLIS UPDATE. # # Arguments: # # initial.x The initial state (a vector) # lpr Function returning the log probability of a state, plus an # arbitrary constant # pf Function returning the random offset (a vector) for a proposal # w The stepsize for proposals, multiplies the offset # initial.lpr The value of lpr(initial.x). May be omitted. # # The value returned is a list containing the following elements: # # next.x The new state # next.lpr The value of lpr(next.x) # rejected TRUE if the proposal was rejected # # The global variable lpr.evals is incremented by the number of calls of # lpr made. metropolis.update <- function (initial.x, lpr, pf, w, initial.lpr=NULL) { # Evaluate the log probability of the initial state, if not already known. if (is.null(initial.lpr)) { initial.lpr <- lpr(initial.x) lpr.evals <<- lpr.evals + 1 } # Propose a candidate state, and evalute its log probability. proposed.x <- initial.x + w*pf() proposed.lpr <- lpr(proposed.x) lpr.evals <<- lpr.evals + 1 # Decide whether to accept or reject the proposed state as the new state. if (runif(1)threshold) { w <- w0 } else { w <- w1 } update <- metropolis.update (states[k,], lpr, pf, w, last.lpr) states[k+1,] <- update$next.x rejected[k] <- update$rejected stepsize[k] <- w rej.in.last.N[1+k%%N] <- update$rejected last.lpr <- update$next.lpr } # Return the states, the rejection flags, and the stepsizes. list (states=states, rejected=rejected, stepsize=stepsize) } # THE SHORT-CUT METROPOLIS METHOD. Simulates a short-cut Metropolis update # consisting of M*L Metropolis updates. # # Arguments: # # initial.x The initial state (a vector of length n) # lpr Function returning the log probability of a state, plus an # arbitrary constant # pf Function returning the random offset (a vector) for a proposal # w The stepsize for proposals, multiplies the offset # L Number of updates in each group # M Number of groups to simulate # min.rej Minimum number of rejections for a good group of L updates # max.rej Maximum number of rejections for a good group of L updates # # The value returned is a list containing the following elements: # # states1 A M*L+1 by n matrix whose rows contain the initial state # and the states after each Metropolis update, including # updates in groups after which a reversal occurred # states2 A M+1 by n matrix whose rows contain the initial state and # the states after each group of L Metropolis updates; the # last row contains the final state from the whole sequence # rejected A vector of length M*L containing the rejection status (T/F) # for each Metropolis update # copied A vector of length M*L containing T/F flags saying whether # each Metropolis update was copied from some earlier state # index A vector of length M*L containing an index for each # Metropolis update that increases from 1 for the first set # of original updates, that decreases from 0 for the second # set of original updates, and is copied for copied updates # # The global variable lpr.evals is incremented by the number of calls of # lpr made. short.cut.metropolis <- function (initial.x, lpr, pf, w, L, M, min.rej=0, max.rej=L-1) { # The function below puts together the list of results that we return. results <- function () { list (states1=states1, states2=states2, rejected=rejected, copied=copied, index=index) } # The function below performs the L Metropolis updates in a group, storing the # results in states1. The last.x and last.lpr variables should contain the # previous state and its log probability; they are updated by this function. # The variable k indexes where in states1 the new states should be stored. # It is incremented in this function. The n.rejected variable is set to the # number of the L updates that were rejections. The rejected, copied, and # index variables are also updated here. The variable ix is incremented by # ixi to keep track of the index for later updates. do.group <- function () { n.rejected <<- 0 for (l in 1:L) { update <- metropolis.update (last.x, lpr, pf, w, last.lpr) states1[k+1,] <<- last.x <<- update$next.x last.lpr <<- update$next.lpr rejected[k] <<- update$rejected copied[k] <<- FALSE index[k] <<- ix <<- ix + ixi k <<- k + 1 n.rejected <<- n.rejected + update$rejected } } # Allocate space for the results. K <- M*L states1 <- matrix(NA,K+1,length(initial.x)) states2 <- matrix(NA,M+1,length(initial.x)) rejected <- rep(NA,K) copied <- rep(NA,K) index <- rep(NA,K) # Store the initial state in the first row of states1 and of states2, and # evaluate its log probability. states1[1,] <- initial.x states2[1,] <- initial.x initial.lpr <- lpr(initial.x) lpr.evals <<- lpr.evals + 1 # Do groups of Metropolis updates starting from the initial state, until we've # done many as were asked for, or until the number of rejections in a group # is outside the limits. States after each Metropolis update are saved # in states1. The final state for each group is saved in states2. Note # that for a group with the number of rejections outside the limits, the # state saved in states2 is not state the last state in states1, but rather # state from the start of the group. The indexes in states1 of the start # (upper0) and end (upper1) of these groups in are recorded for later use # in copying states. last.x <- initial.x last.lpr <- initial.lpr k <- 1 ix <- 0 ixi <- 1 upper0 <- k while (k<=K) { states2[2+(k-1)/L,] <- last.x do.group() if (n.rejectedmax.rej) { upper1 <- k break } else { states2[1+(k-1)/L,] <- states1[k,] } } if (k>K) return (results()) # Return if we've done all M groups. # Copy already-computed states as "new" states as we move backwards through # the groups previously simulated. Don't copy the last group causing the # reversal, for which the number of rejections was outside the limits. j <- upper1 - L - 1 while (k<=K && j>=upper0) { states1[(k+1):(k+L),] <- states1[j:(j-L+1),] rejected[k:(k+L-1)] <- rejected[j:(j-L+1)] copied[k:(k+L-1)] <- TRUE index[k:(k+L-1)] <- index[j:(j-L+1)] states2[1+(k+L-1)/L,] <- states2[1+(j-L)/L,] k <- k + L j <- j - L } if (k>K) return (results()) # Return if we've done all M groups. # Restore the initial state, then do more groups of Metropolis updates, # until we've done as many as were asked for, or the number of rejections in # a group is outside the limits. Record the indexes of the start (lower0) # and end (lower1) of these groups in states1. last.x <- initial.x last.lpr <- initial.lpr ix <- 1 ixi <- -1 lower0 <- k while (k<=K) { states2[2+(k-1)/L,] <- last.x do.group() if (n.rejectedmax.rej) { lower1 <- k break } else { states2[1+(k-1)/L,] <- states1[k,] } } if (k>K) return (results()) # Return if we've done all M groups. # Copy already-computed states as "new" states, going back and forth over # the "lower" groups and the "upper" groups. repeat { # Copy the lower states backwards, excluding the group causing the reversal. j <- lower1 - L - 1 while (k<=K && j>=lower0) { states1[(k+1):(k+L),] <- states1[j:(j-L+1),] rejected[k:(k+L-1)] <- rejected[j:(j-L+1)] index[k:(k+L-1)] <- index[j:(j-L+1)] copied[k:(k+L-1)] <- TRUE states2[1+(k+L-1)/L,] <- states2[1+(j-L)/L,] k <- k + L j <- j - L } if (k>K) return (results()) # Return if we've done all M groups. # Copy the upper states forwards, including the group causing the reversal. j <- upper0+1 while (k<=K && j<=upper1) { states1[(k+1):(k+L),] <- states1[j:(j+L-1),] rejected[k:(k+L-1)] <- rejected[(j-1):(j+L-2)] index[k:(k+L-1)] <- index[(j-1):(j+L-2)] copied[k:(k+L-1)] <- TRUE states2[1+(k+L-1)/L,] <- states2[1+(j+L-2)/L,] k <- k + L j <- j + L } if (k>K) return (results()) # Return if we've done all M groups. # Copy the upper states backwards, excluding the group causing the reversal. j <- upper1 - L - 1 while (k<=K && j>=upper0) { states1[(k+1):(k+L),] <- states1[j:(j-L+1),] rejected[k:(k+L-1)] <- rejected[j:(j-L+1)] index[k:(k+L-1)] <- index[j:(j-L+1)] copied[k:(k+L-1)] <- TRUE states2[1+(k+L-1)/L,] <- states2[1+(j-L)/L,] k <- k + L j <- j - L } if (k>K) return (results()) # Return if we've done all M groups. # Copy the lower states forwards, including the group causing the reversal. j <- lower0 + 1 while (k<=K && j<=lower1) { states1[(k+1):(k+L),] <- states1[j:(j+L-1),] rejected[k:(k+L-1)] <- rejected[(j-1):(j+L-2)] index[k:(k+L-1)] <- index[(j-1):(j+L-2)] copied[k:(k+L-1)] <- TRUE states2[1+(k+L-1)/L,] <- states2[1+(j+L-2)/L,] k <- k + L j <- j + L } if (k>K) return (results()) # Return if we've done all M groups. } } # PERFORM SEVERAL SHORT-CUT METROPOLIS SEQUENCES. # # Arguments: # # initial.x The initial state (a vector of length n) # lpr Function returning the log probability of a state, plus an # arbitrary constant # pf Function returning the random offset (a vector) for a proposal # w.vec Vector of stepsizes for proposals # L.vec Vector of numbers of updates in each group # M.vec Vector of numbers of groups # min.rej.vec Vector of minimum numbers of rejections for a good group # max.rej.vec Vector of maximum numbers of rejections for a good group # # Returns the concatenation of the results of short.cut.metropolis with w, L, # M, min.rej, and max.rej taking on the values in w.vec, L.vec, M.vec, # min.rej.vec, and max.rej.vec (which are extended by repetition to the length # of the longest). In detail, the value returned is a list with the following # elements: # # states1 A matrix whose rows contain the initial state and the states # after each Metropolis update, including updates in groups # after which a reversal occurred # states2 A matrix whose rows contain the initial state and the states # after each group of L Metropolis updates; the last row # contains the final state after all short-cut sequences # rejected A vector containing the rejection status (T/F) of each update # copied A vector containing T/F flags saying whether each update # was copied from some earlier state # index A vector containing an index for each update that increases # from 1 for the first set of original updates, that decreases # from 0 for the second set of original updates, and is copied # for copied updates # stepsize A vector containing the stepsize used for each update # # The global variable lpr.evals is incremented by the number of calls of # lpr made. short.cut.sequences <- function (initial.x, lpr, pf, w.vec, L.vec, M.vec, min.rej.vec=0, max.rej.vec=L.vec-1) { N <- max (length(w.vec), length(L.vec), length(M.vec), length(min.rej.vec), length(max.rej.vec)) w.vec <- rep(w.vec,length=N) L.vec <- rep(L.vec,length=N) M.vec <- rep(M.vec,length=N) min.rej.vec <- rep(min.rej.vec,length=N) max.rej.vec <- rep(max.rej.vec,length=N) states1 <- matrix(NA,sum(L.vec*M.vec)+1,length(initial.x)) states2 <- matrix(NA,sum(M.vec)+1,length(initial.x)) rejected <- c() copied <- c() index <- c() stepsize <- c() states1[1,] <- initial.x states2[1,] <- initial.x k <- 1 m <- 1 for (i in 1:N) { res <- short.cut.metropolis (states2[m,], lpr, pf, w.vec[i], L.vec[i], M.vec[i], min.rej.vec[i], max.rej.vec[i]) states1[k:(k+L.vec[i]*M.vec[i]),] <- res$states1 states2[m:(m+M.vec[i]),] <- res$states2 rejected <- c(rejected,res$rejected) copied <- c(copied,res$copied) index <- c(index,res$index) stepsize <- c(stepsize,rep(w.vec[i],L.vec[i]*M.vec[i])) k <- k + L.vec[i]*M.vec[i] m <- m + M.vec[i] } list (states1=states1, states2=states2, rejected=rejected, copied=copied, index=index, stepsize=stepsize) } # PERFORM SEVERAL STANDARD METROPOLIS SEQUENCES. # # Arguments: # # initial.x The initial state (a vector of length n) # lpr Function returning the log probability of a state, plus an # arbitrary constant # pf Function returning the random offset (a vector) for a proposal # w.vec Vector of stepsizes for proposals # K.vec Vector of number of updates for each stepsize # # Returns the concatenation of the results of standard.metropolis with w and K # taking on the values in w.vec and K.vec (which are extended by repetition to # the length of the longest). In detail, the value returned is a list with the # following elements: # # states A matrix whose rows contain the initial state and the states # after each Metropolis update. # rejected A vector containing the rejection status (T/F) of each update # stepsize A vector containing the stepsize used for each update # # The global variable lpr.evals is incremented by the number of calls of # lpr made. standard.sequences <- function (initial.x, lpr, pf, w.vec, K.vec) { N <- max (length(w.vec), length(K.vec)) w.vec <- rep(w.vec,length=N) K.vec <- rep(K.vec,length=N) states <- matrix(NA,sum(K.vec)+1,length(initial.x)) rejected <- c() stepsize <- c() states[1,] <- initial.x k <- 1 for (i in 1:N) { res <- standard.metropolis (states[k,], lpr, pf, w.vec[i], K.vec[i]) states[k:(k+K.vec[i]),] <- res$states rejected <- c(rejected,res$rejected) stepsize <- c(stepsize,rep(w.vec[i],K.vec[i])) k <- k + K.vec[i] } list (states=states, rejected=rejected, stepsize=stepsize) } # PERFORM SEVERAL SHORT-CUT METROPOLIS UPDATES, RETAINING ONLY FINAL STATES. # # Arguments: # # initial.x The initial state (a vector of length n) # lpr Function returning the log probability of a state, plus an # arbitrary constant # pf Function returning the random offset (a vector) for a proposal # w.vec Vector of stepsizes for proposals # L.vec Vector of numbers of updates in each group # M.vec Vector of numbers of groups # min.rej.vec Vector of minimum numbers of rejections for a good group # max.rej.vec Vector of maximum numbers of rejections for a good group # # Forms a matrix whose rows contain the initial state and the final states # from short.cut.metropolis with w, L, M, min.rej, and max.rej taking on the # values in w.vec, L.vec, M.vec, min.rej.vec, and max,rej.vec (which are # extended by repetition to the length of the longest). The states after # each update and after each group of updates (except for the final state) # are discarded. In detail, the value returned is a list with the following # elements: # # states A matrix whose rows contain the initial state and the final # states of each short-cut sequence # rejected A vector containing the fraction of rejections for each # short-cut sequence # copied A vector containing the fraction of copied states for each # short-cut sequence # stepsize A vector containing the stepsize used for each short-cut # sequence # repeated A vector containing the number of Metropolis updates in # each short-cut sequence short.cut.final <- function (initial.x, lpr, pf, w.vec, L.vec, M.vec, min.rej.vec=0, max.rej.vec=L.vec-1) { N <- max (length(w.vec), length(L.vec), length(M.vec), length(min.rej.vec), length(max.rej.vec)) w.vec <- rep(w.vec,length=N) L.vec <- rep(L.vec,length=N) M.vec <- rep(M.vec,length=N) min.rej.vec <- rep(min.rej.vec,length=N) max.rej.vec <- rep(max.rej.vec,length=N) states <- matrix(NA,N+1,length(initial.x)) rejected <- c() copied <- c() states[1,] <- initial.x k <- 1 for (i in 1:N) { res <- short.cut.metropolis (states[k,], lpr, pf, w.vec[i], L.vec[i], M.vec[i], min.rej.vec[i], max.rej.vec[i]) states[k+1,] <- res$states2[nrow(res$states2),] rejected <- c(rejected,mean(as.numeric(res$rejected))) copied <- c(copied,mean(as.numeric(res$copied))) k <- k + 1 } list (states=states, rejected=rejected, copied=copied, stepsize=w.vec, repeated=L.vec*M.vec) } # PERFORM SEVERAL STANDARD METROPOLIS SEQUENCES, RETAINING ONLY FINAL STATES. # # Arguments: # # initial.x The initial state (a vector of length n) # lpr Function returning the log probability of a state, plus an # arbitrary constant # pf Function returning the random offset (a vector) for a proposal # w.vec Vector of stepsizes for proposals # K.vec Vector of numbers of updates to do with each stepsize # # Forms a matrix whose rows contain the initial state and the final states # from performing sequences of standard Metropolis updates with w and K taking # on the values in w.vec and K.vec (which are extended by repetition to the # length of the longest). In detail, the value returned is a list with the # following elements: # # states A matrix whose rows contain the initial state and the states # after each complete sequence of K updates # rejected A vector containing the fraction of rejections for each # sequence # repeated A vector containing the number of Metropolis updates in # each sequence # stepsize A vector containing the stepsize used for each sequence standard.final <- function (initial.x, lpr, pf, w.vec, K.vec) { N <- max(length(w.vec),length(K.vec)) w.vec <- rep(w.vec,length=N) K.vec <- rep(K.vec,length=N) states <- matrix(NA,N+1,length(initial.x)) rejected <- c() states[1,] <- initial.x k <- 1 for (i in 1:N) { res <- standard.metropolis (states[k,], lpr, pf, w.vec[i], K.vec[i]) states[k+1,] <- res$states[nrow(res$states),] rejected <- c(rejected,mean(as.numeric(res$rejected))) k <- k + 1 } list (states=states, rejected=rejected, stepsize=w.vec, repeated=K.vec) }