In this lab we’ll simulate some “random walks” as a way to generate samples from a target distribution. We’ll introduce the Metropolis Algorithm and see how it can be used to “shape” the samples that our random walk gives us.
As an example goal, suppose we want to generate samples from a random process that returns integer values between 1 and 6, but unlike a fair die, the values are not equally likely. Instead, they should occur with probabilities given by the PMF:
\[ p(\theta) = \theta / 21 \qquad \theta = 1, 2, 3, 4, 5, 6\] You should verify that these probabilities sum to 1.
A graph of this PMF looks like so:
<- tibble(theta = 1:6) %>%
PMFTable mutate(p_theta = theta / 21)
%>%
PMFTable ggplot(aes(x = factor(theta), y = p_theta)) +
geom_bar(stat = "identity", fill = "dodgerblue") +
scale_x_discrete(
name = expression(theta),
breaks = 1:6
+
) scale_y_continuous(
name = expression(p(theta)),
limits = c(0,0.5),
breaks = seq(0, 0.5, by = 0.02),
minor_breaks = NULL
)
This is a simple enough distribution that it’s easy to sample from directly. For most real posterior distributions we won’t be able to do this direct sampling, but I think it’s useful to do a comparison using a case where we can sample directly so we can see how the estimation error for expected value relates to the number of samples we have.
Here’s code that generates 1000 samples from this distribution.
set.seed(8675309)
<- tibble(
DirectSamples s = 1:1000,
theta = sample(1:6, size = 1000, replace = TRUE, prob = (1:6)/21))
Here’s a plot of the results by “iteration”, with successive samples connected by lines (this is just for comparison to what happens with sequential samples)
<- DirectSamples %>%
path_plot gf_line(theta ~ s, color = "tomato") +
scale_x_continuous(
name = "Iteration",
breaks = seq(0, 1000, by = 50)
+
) scale_y_continuous(
name = expression(theta^(s)),
breaks = 1:6
) path_plot
and here’s a plot of how often we got each value. The actual PMF values are shown by blue dots for reference.
<- DirectSamples %>%
frequency_plot gf_props(~factor(theta), fill = "tomato") +
geom_point(aes(x = theta, y = p_theta), data = PMFTable, color = "dodgerblue") +
scale_x_discrete(
name = expression(theta),
breaks = 1:6
+
) scale_y_continuous(
name = expression(p(theta)),
limits = c(0,0.5),
breaks = seq(0, 0.5, by = 0.02),
minor_breaks = NULL
) frequency_plot
Recall that we could define a random walk around the 1-6 “circle” by flipping two coins, moving “counterclockwise” if we get two tails, staying put if we get one head and one tails, and moving “clockwise” if we get two heads.
After enough iterations of this, we saw that the process’s “memory” of where it started fades away, and we wind up equally likely to be at any one of the six values.
Here’s a simulation of this process in code:
set.seed(29747)
<- tibble(s = NULL, theta = NULL)
RandomWalk ## Initialize our position at iteration s = 0 to theta = 1
<- 1
theta for(iter in 1:1000) {
## Generate the step with 1/4 chance of going left, 1/2 chance of staying
## and 1/4 chance of going right
<- sample(c(-1, 0, 1), size = 1, prob = c(0.25, 0.50, 0.25))
step ## add the step to the last theta
<- theta + step
theta ## wrap around
if(theta == 0) theta <- 6
if(theta == 7) theta <- 1
## Record the result
"s"] <- iter
RandomWalk[iter, "theta"] <- theta
RandomWalk[iter, }
Here’s a plot of the results by iteration:
%>%
RandomWalk gf_line(theta ~ s, color = "tomato") +
scale_x_continuous(
name = "Iteration",
breaks = seq(0, 1000, by = 50)
+
) scale_y_continuous(
name = expression(theta^(s)),
breaks = 1:6
)
and here’s a plot of how often we end up getting each value:
%>%
RandomWalk gf_props(~factor(theta), fill = "tomato") +
geom_point(aes(x = theta, y = 1/6), data = PMFTable, color = "dodgerblue") +
scale_x_discrete(
name = expression(theta),
breaks = 1:6
+
) scale_y_continuous(
name = expression(p(theta)),
limits = c(0,0.5),
breaks = seq(0, 0.5, by = 0.02),
minor_breaks = NULL
)
For comparison, here’s a set of 1000 direct independent samples from the uniform distribution.
<- tibble(
DirectSamples_Uniform s = 1:1000,
theta = sample(1:6, size = 1000, replace = TRUE, prob = rep(1/6, 6)))
%>%
DirectSamples_Uniform gf_line(theta ~ s, color = "tomato") +
scale_x_continuous(
name = "Iteration",
breaks = seq(0, 1000, by = 50)
+
) scale_y_continuous(
name = expression(theta^(s)),
breaks = 1:6
)
%>%
DirectSamples_Uniform gf_props(~factor(theta), fill = "tomato") +
geom_point(aes(x = theta, y = 1/6), data = PMFTable, color = "dodgerblue") +
scale_x_discrete(
name = expression(theta),
breaks = 1:6
+
) scale_y_continuous(
name = expression(p(theta)),
limits = c(0,0.5),
breaks = seq(0, 0.5, by = 0.02),
minor_breaks = NULL
)
This is helpful if creating a uniform distribution is our goal, but can we modify this process so that we wind up with probabilities that match our target distribution?
The key will be that we want to introduce asymmetric “friction” to our random walk to cause it to “stick” at values that we want to visit more often and “resist” moving to values that we want to visit less often.
We can achieve this by introducing a probabilistic “rejection” step at each iteration: rather than automatically moving according to our coin flips, we introduce a chance of staying where we are if the place we are now should be more probable than the place we’re trying to go.
In a nutshell, what we’ll do is to break up each iteration into a few parts. Assuming our current location is \(\theta^{(s-1)}\), we will:
Here’s some code that modifies our previous random walk to include this “accept/reject” step.
set.seed(29747)
## Initialize our position at iteration s = 0 to theta = 1
<- 1
theta <- tibble(s = 0, theta = theta)
MetropolisWalk ## Define a function that gives the target PMF
<- function(theta) {return(theta / 21)}
p_theta for(iter in 1:1000) {
## Generate the step with 1/4 chance of going left, 1/2 chance of staying
## and 1/4 chance of going right
<- sample(c(-1, 0, 1), size = 1, prob = c(0.25, 0.50, 0.25))
step ## add the step to the last theta
<- theta + step
theta_star ## wrap around
if(theta_star == 0) theta_star <- 6
if(theta_star == 7) theta_star <- 1
## calculate the acceptance probability
<- min(1, p_theta(theta_star) / p_theta(theta))
alpha ## determine whether to accept the proposal with the corresponding probability
if(alpha < 1) {
<- sample(c(0,1), prob = c(1 - alpha, alpha))
accept else {
} <- 1
accept
}## set the new theta accordingly
if(accept == 1) theta <- theta_star
else theta <- theta
## Record the result
"s"] <- iter
MetropolisWalk[iter, "theta"] <- theta
MetropolisWalk[iter, }
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used
## Warning in if (accept == 1) theta <- theta_star else theta <- theta: the
## condition has length > 1 and only the first element will be used