# SHORT-MET-ABBREV.R - Abbreviated version of short-met.r.
# 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).
#
# 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
metropolis.update <- function (initial.x, lpr, pf, w, initial.lpr)
{
# Propose a candidate state, and evalute its log probability.
proposed.x <- initial.x + w*pf()
proposed.lpr <- lpr(proposed.x)
# Decide whether to accept or reject the proposed state as the new state.
if (runif(1)max.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),]
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
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),]
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),]
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),]
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),]
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.
}
}