prob <- 0.1
len_z <- 1000
z <- rbinom(n = len_z, 1, prob)
n <- 1000

z_mat <- matrix(nrow = n, ncol = len_z)

for(i in 1:n){
  for(s in 1:len_z){
    oldz <- z[s]
    newz <- rbinom(1, 1, 0.5)

    alpha <- min(1, (prob^newz * (1 - prob)^(1 - newz))/(prob^oldz * (1 - prob)^(1 - oldz)))

    if(rbinom(1, 1, alpha) == 1){
      z[s] <- newz
    }
  }
  z_mat[i,] <- z
}

x = apply(z_mat, 1, mean)
mean(x)