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