Figure 7.3

[ASKurz]

In Section 7.2.3, General properties of a random walk, Kruschke covered Figure 7.3. From the beginning of the section, we read:

The trajectory shown in Figure 7.2 is just one possible sequence of positions when the movement heuristic is applied. At each time step, the direction of the proposed move is random, and if the relative probability of the proposed position is less than that of the current position, then acceptance of the proposed move is also random. Because of the randomness, if the process were started over again, then the specific trajectory would almost certainly be different. Regardless of the specific trajectory, in the long run the relative frequency of visits mimics the target distribution.

Figure 7.3 shows the probability of being in each position as a function of time. At time t = 1, the politician starts at θ = 4. This starting position is indicated in the upper- left panel of Figure 7.3, labeled t = 1, by the fact that there is 100% probability of being at θ = 4.

We want to derive the probability of ending up in each position at the next time step. To determine the probabilities of positions for time t = 2, consider the possibilities from the movement process. The process starts with the flip of a fair coin to decide which direction to propose moving. (p. 149)

We see the figure on the next page:

At present, I’m not sure how to pull off the simulation necessary to generate the figure. If you have the random-walk simulation chops, please share. If at all possible tidyverse-oriented workflows are preferred.

[cmoten]

Dr. Kruschke made a blog entry at this link. In essence, you have to create a transition matrix that is composed of a matrix of proposal and acceptance probabilities. You will then create a vector of position probabilities at each time period. This forms the Markov chain that forms the long-run transition probabilities which will reach the target probabilities as the time approaches infinity.

I recreated the plot of this post by using this code:

library(ggplot2)
library(gridExtra)

nslots <- 7
p_target <- 1:7
p_target <- p_target/sum(p_target)

#Construct the transition matrix

proposal_matrix <- matrix(0, nrow = nslots, ncol = nslots)

for(from_idx in 1:nslots){
  for(to_idx in 1:nslots) {
    if(to_idx == from_idx -1) {proposal_matrix[from_idx, to_idx] <- 0.5}
    if(to_idx == from_idx + 1) {proposal_matrix[from_idx, to_idx] <- 0.5}
  }
}

#Construct the acceptance matrix
acceptance_matrix <- matrix(0, nrow = nslots, ncol = nslots)

for(from_idx in 1:nslots) {
  for(to_idx in 1:nslots) {
    acceptance_matrix[from_idx, to_idx] <- min(p_target[to_idx]/p_target[from_idx], 1)
  }
}

#Compute the matrix of move probabilities
move_matrix <- proposal_matrix * acceptance_matrix

# Compute the transition matrix, including the probability of staying in place:
transition_matrix <- move_matrix
for ( diag_idx in 1:nslots ) {
  transition_matrix[diag_idx, diag_idx] = 1.0 - sum(move_matrix[diag_idx,])
}
show(transition_matrix)
#>      [,1]      [,2]      [,3]  [,4]      [,5]       [,6]      [,7]
#> [1,] 0.50 0.5000000 0.0000000 0.000 0.0000000 0.00000000 0.0000000
#> [2,] 0.25 0.2500000 0.5000000 0.000 0.0000000 0.00000000 0.0000000
#> [3,] 0.00 0.3333333 0.1666667 0.500 0.0000000 0.00000000 0.0000000
#> [4,] 0.00 0.0000000 0.3750000 0.125 0.5000000 0.00000000 0.0000000
#> [5,] 0.00 0.0000000 0.0000000 0.400 0.1000000 0.50000000 0.0000000
#> [6,] 0.00 0.0000000 0.0000000 0.000 0.4166667 0.08333333 0.5000000
#> [7,] 0.00 0.0000000 0.0000000 0.000 0.0000000 0.42857143 0.5714286

# Specify starting position vector:
position_vec = rep(0,nslots)
position_vec[round(nslots/2)] = 1.0

p <- list()
data <- as.data.frame(cbind(position = 1:nslots, prob = position_vec))

#Loop through the requisite time indexes
#Create a plot
#Update the data and transition vector
for(time_idx in 1:15) {
  #browser()
  p[[time_idx]] <- ggplot(data, aes(x = position, y = prob)) +
    geom_segment(aes(x = position, y = 0, xend = position, yend = prob)) +
    scale_x_continuous(expression(theta), breaks = 1:7) +
    ylab(expression(paste("p(", theta, ")"))) +
    ggtitle(bquote(t==.(time_idx))) +
    theme_minimal()

  #Update the position vec
  position_vec <- position_vec %*% transition_matrix

  data <- NULL #In order to reset the data for the next plot
  data <- as.data.frame(cbind(position = 1:nslots, prob = t(position_vec)))
  names(data) <- c("position", "prob")
}

#Plot the target distribution
data <- NULL
data <- data.frame(position = 1:nslots, prob = p_target)
p[[time_idx + 1]] <- ggplot(data, aes(x = position, y = prob)) +
  geom_segment(aes(x = position, y = 0, xend = position, yend = prob)) +
  scale_x_continuous(expression(theta), breaks = 1:7) +
  ylab(expression(paste("p(", theta, ")"))) +
  ggtitle("target") +
  theme_minimal()

#Layout the final plot
do.call("grid.arrange", c(p, nrow = 4, as.table = FALSE))

^{Created on 2020-04-03 by the reprex package (v0.3.0)}

[ASKurz]

@cmoten, I don't have time to look through this in detail, today. But on first glance, this is excellent. Thank you so much!

[cmoten]

@ASKurz no problem and thank you for providing this book. I am interested in learning more about the brms and tidybayes packages and your book provides a good start to learn them.

[ASKurz]

Finally had a chance to walk through this. It works! I’ll update this section and credit you for your answer in the next release of the book.