```

References

Load packages

library(magrittr)
library(assertthat)
library(ggplot2)
library(tidyr)

Constant \(k\) is taken to be one.

position <- function(U1) {
    cut(U1, breaks = c(0,0.25, 0.5, 0.75, 1.0), labels = FALSE, include.lowest = TRUE)
}

prob_plus <- function(k, neighbors) {
    p_plus  <- prod(exp(k * +1 * neighbors))
    p_minus <- prod(exp(k * -1 * neighbors))
    ## Normalize
    p_plus / (p_plus + p_minus)
}

update_v <- function(v, k, U1, U2) {
    assert_that(U1 >= 0)
    assert_that(U1 <= 1)
    assert_that(U2 >= 0)
    assert_that(U2 <= 1)
    assert_that(length(v) == 4)
    assert_that(length(k) == 1)
    assert_that(length(U1) == 1)
    assert_that(length(U2) == 1)

    ## Choose target position
    pos <- position(U1)

    if (pos == 1) {
        neighbors <- c(v[4],v[2])
    } else if (pos == 2) {
        neighbors <- c(v[1],v[3])
    } else if (pos == 3) {
        neighbors <- c(v[2],v[4])
    } else if (pos == 4) {
        neighbors <- c(v[3],v[1])
    }

    ## Update target position based on neighbors (using Gibbs sampling)
    v[pos] <- ifelse(U2 < prob_plus(k = k, neighbors = neighbors),
                     +1, -1)

    v
}

v_max <- c(+1,+1,+1,+1)
v_min <- c(-1,-1,-1,-1)

update_v(v_max, k = 1, 0.5, 0.5)
## [1] 1 1 1 1
n <- 100
U1_seq <- runif(n = n)
U2_seq <- runif(n = n)

run_chain <- function(init, k, U1_seq, U2_seq) {
    assert_that(length(init) == 4)
    assert_that(length(U1_seq) == length(U2_seq))

    ## Out object
    out <- matrix(NA,
                  nrow = length(U1_seq),
                  ncol = length(init))
    out[1,] <- init

    for (i in seq_along(U1_seq)[-1]) {
        out[i,] <- update_v(out[i-1,], k, U1_seq[i], U2_seq[i])
    }

    out
}

status <- function(output){
  assert_that(is.matrix(output) == TRUE)
  output[output == -1] <- 0
  out <- rep(NA, nrow(output))
  out <- apply(output, 1, function(row){2^3 * row[1] + 2^2 * row[2] + 2^1 * row[3] + row[4]})
  out
}

Try different \(k\)

With a smaller \(k\), convergence seems faster.

data1 <- data.frame(i = seq_along(U1_seq),
                    max_chain_k0.01 = rowSums(run_chain(v_max, 0.01, U1_seq, U2_seq)),
                    min_chain_k0.01 = rowSums(run_chain(v_min, 0.01, U1_seq, U2_seq)),
                    max_chain_k0.5 = rowSums(run_chain(v_max, 0.5, U1_seq, U2_seq)),
                    min_chain_k0.5 = rowSums(run_chain(v_min, 0.5, U1_seq, U2_seq)),
                    max_chain_k1 = rowSums(run_chain(v_max, 1, U1_seq, U2_seq)),
                    min_chain_k1 = rowSums(run_chain(v_min, 1, U1_seq, U2_seq)),
                    max_chain_k2 = rowSums(run_chain(v_max, 2, U1_seq, U2_seq)),
                    min_chain_k2 = rowSums(run_chain(v_min, 2, U1_seq, U2_seq)))

data1_long <- gather(data1, key = chain, value = value, -i)
data1_long$k <- as.numeric(gsub(".*_chain_k", "", data1_long$chain))

ggplot(data = data1_long,
       mapping = aes(x = i, y = value, group = chain, color = chain)) +
    geom_line(alpha = 1/3, size = 2) +
    facet_grid(k ~ .) +
    theme_bw() + theme(legend.key = element_blank())

Draw the actual path over status 0 - 15 using binary expansion in order to verify MAX and MIN indeed converge into the same distribution.

data1_status <- data.frame(i = seq_along(U1_seq),
                    max_chain_k0.01 = status(run_chain(v_max, 0.01, U1_seq, U2_seq)),
                    min_chain_k0.01 = status(run_chain(v_min, 0.01, U1_seq, U2_seq)),
                    max_chain_k0.5 = status(run_chain(v_max, 0.5, U1_seq, U2_seq)),
                    min_chain_k0.5 = status(run_chain(v_min, 0.5, U1_seq, U2_seq)),
                    max_chain_k1 = status(run_chain(v_max, 1, U1_seq, U2_seq)),
                    min_chain_k1 = status(run_chain(v_min, 1, U1_seq, U2_seq)),
                    max_chain_k2 = status(run_chain(v_max, 2, U1_seq, U2_seq)),
                    min_chain_k2 = status(run_chain(v_min, 2, U1_seq, U2_seq)))

data1_long_status <- gather(data1_status, key = chain, value = value, -i)
data1_long_status$k <- as.numeric(gsub(".*_chain_k", "", data1_long$chain))

ggplot(data = data1_long_status,
       mapping = aes(x = i, y = value, group = chain, color = chain)) +
    geom_line(alpha = 1/3, size = 2) +
    facet_grid(k ~ .) +
    theme_bw() + theme(legend.key = element_blank())

Try starting at different index

\(k = 0.01\) for faster converence. No matter when you start, convergence into the same distribution occurs.

## Try different starting index
get_sum <- function(init, k, start_i, ran_data, status = FALSE) {
    assert_that(start_i < nrow(ran_data))
    n <- nrow(ran_data)
    assert_that(start_i >= 1)
    if (status == FALSE){
      c(rep(NA, start_i - 1),
        rowSums(run_chain(init,
                        k,
                        ran_data[seq(start_i, n),"U1_seq"],
                        ran_data[seq(start_i, n),"U2_seq"])))
    }else{
      c(rep(NA, start_i - 1),
        status(run_chain(init,
                        k,
                        ran_data[seq(start_i, n),"U1_seq"],
                        ran_data[seq(start_i, n),"U2_seq"])))
    }
}

n <- 100

ran_data <- data.frame(i = seq_len(n),
                       U1_seq = runif(n = n),
                       U2_seq = runif(n = n))

data2 <- data.frame(i = seq_len(nrow(ran_data)),
                    max_chain_i1  = get_sum(v_max, k = 0.01, start_i = 1, ran_data = ran_data),
                    min_chain_i1  = get_sum(v_min, k = 0.01, start_i = 1, ran_data = ran_data),
                    max_chain_i20 = get_sum(v_max, k = 0.01, start_i = 20, ran_data = ran_data),
                    min_chain_i20 = get_sum(v_min, k = 0.01, start_i = 20, ran_data = ran_data),
                    max_chain_i30 = get_sum(v_max, k = 0.01, start_i = 30, ran_data = ran_data),
                    min_chain_i30 = get_sum(v_min, k = 0.01, start_i = 30, ran_data = ran_data),
                    max_chain_i40 = get_sum(v_max, k = 0.01, start_i = 40, ran_data = ran_data),
                    min_chain_i40 = get_sum(v_min, k = 0.01, start_i = 40, ran_data = ran_data))


data2_long <- gather(data2, key = chain, value = value, -i)
data2_long$start_i <- as.numeric(gsub(".*_chain_i", "", data2_long$chain))

ggplot(data = data2_long,
       mapping = aes(x = i, y = value, group = chain, color = chain)) +
    geom_line(alpha = 1/3, size = 2) +
    theme_bw() + theme(legend.key = element_blank())
## Warning: Removed 174 rows containing missing values (geom_path).

data2_status <- data.frame(i = seq_len(nrow(ran_data)),
                    max_chain_i1  = get_sum(v_max, k = 0.01, start_i = 1, ran_data = ran_data, status = TRUE),
                    min_chain_i1  = get_sum(v_min, k = 0.01, start_i = 1, ran_data = ran_data, status = TRUE),
                    max_chain_i20 = get_sum(v_max, k = 0.01, start_i = 20, ran_data = ran_data, status = TRUE),
                    min_chain_i20 = get_sum(v_min, k = 0.01, start_i = 20, ran_data = ran_data, status = TRUE),
                    max_chain_i30 = get_sum(v_max, k = 0.01, start_i = 30, ran_data = ran_data, status = TRUE),
                    min_chain_i30 = get_sum(v_min, k = 0.01, start_i = 30, ran_data = ran_data, status = TRUE),
                    max_chain_i40 = get_sum(v_max, k = 0.01, start_i = 40, ran_data = ran_data, status = TRUE),
                    min_chain_i40 = get_sum(v_min, k = 0.01, start_i = 40, ran_data = ran_data, status = TRUE))


data2_long_status <- gather(data2_status, key = chain, value = value, -i)
data2_long_status$start_i <- as.numeric(gsub(".*_chain_i", "", data2_long$chain))

ggplot(data = data2_long_status,
       mapping = aes(x = i, y = value, group = chain, color = chain)) +
    geom_line(alpha = 1/3, size = 2) +
    theme_bw() + theme(legend.key = element_blank())
## Warning: Removed 174 rows containing missing values (geom_path).

Another toy example

We also consider a toy example showing that the converged distribution (coupling from the past) is indeed the stationary distribution \(\pi\) of the Markov Chain, which we did not verify directly in the above example. Here we consider the transition matrix to be \[P = \left(\begin{array}{cc}0.5 & 0.5 \\ 0.4 & 0.6\end{array}\right).\]

We can solve that the stationary distribution that satisfies \(\pi P = \pi\) is \[\pi = \left(\frac{4}{9}, \frac{5}{9}\right).\]

update_v <-function(v, U){
  if (v == 0 & U < 0.5){
  v <- 0
  }else if (v == 0 & U >= 0.5){
  v <- 1
  }else if (v == 1 & U < 0.4){
  v <- 0
  }else{v <- 1}
  v
}

couple <- function(min = 0, max = 1, U_seq) {
  chain1 <- min
  chain2 <- max
  for (i in 1:length(U_seq)) {
    chain1 <- update_v(chain1, U_seq[i])
    chain2 <- update_v(chain2, U_seq[i])
  }
  if (chain1 != chain2){return(NA)}else{return(chain1)}
}

run_chain <- function(min = 0, max = 1) {
  U_seq <- runif(1,0,1)
  result <- NA
  while (is.na(result)) {
    result <- couple(min = min, max = max, U_seq)
    ## Trick: Shuffle the order of U(0,1)
    U_seq <- c(runif(1,0,1), U_seq)  
  }
  return(c(result, length(U_seq - 1)))
}

n <- 10000
output <- matrix(NA, nrow = n, ncol = 2)
for (k in 1:n){
  output[k, ]<- run_chain()
}

data3 <- data.frame(value = output[,1])
mean(data3$value)
## [1] 0.5521
## Compare to the theoretical value.
5/9
## [1] 0.5555556
chisq.test(table(data3$value), p = c(4/9, 5/9))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(data3$value)
## X-squared = 0.48361, df = 1, p-value = 0.4868

Thus, once the MAX and MIN chains coupled, it is indeed converges to the stationary distribution \(\pi\).