加速一维随机游走模拟：提升马尔可夫链仿真效率的方法咨询

问题背景 #

一个变量以0.5的概率+1，0.5的概率-1移动（一维随机游走马尔可夫链），起始位置为5。

场景1 #

求该变量首次到达位置0的期望时间？

场景2 #

求该变量在至少到达过位置10后，首次到达位置0的期望时间？

尝试的模拟代码（运行耗时过长已中断） #

# the simulations can take a long time to run, I interrupted them

library(ggplot2)
library(gridExtra)

n_sims <- 100
times_to_end_0 <- numeric(n_sims)
times_to_end_0_after_10 <- numeric(n_sims)
paths_0 <- vector("list", n_sims)
paths_0_after_10 <- vector("list", n_sims)

for (i in 1:n_sims) {
    print(paste("Running simulation", i, "for situation 1..."))
    y <- 5
    time <- 0
    path_0 <- c(y)

    while(y > 0) {
        step <- sample(c(-1, 1), 1, prob = c(0.5, 0.5))
        y <- y + step
        path_0 <- c(path_0, y)
        time <- time + 1

        if (y == 0) {
            times_to_end_0[i] <- time
            paths_0[[i]] <- data.frame(time = 1:length(path_0), y = path_0, sim = i)
            break
        }
    }

    print(paste("Running simulation", i, "for situation 2..."))
    y <- 5
    time <- 0
    reached_10 <- FALSE
    path_0_after_10 <- c(y)

    while(y > 0 || !reached_10) {
        step <- sample(c(-1, 1), 1, prob = c(0.5, 0.5))
        y <- y + step
        path_0_after_10 <- c(path_0_after_10, y)
        time <- time + 1

        if (y == 10) {
            reached_10 <- TRUE
        }

        if (y == 0 && reached_10) {
            times_to_end_0_after_10[i] <- time
            paths_0_after_10[[i]] <- data.frame(time = 1:length(path_0_after_10), y = path_0_after_10, sim = i)
            break
        }
    }
}

df1 <- data.frame(time = times_to_end_0)
df2 <- data.frame(time = times_to_end_0_after_10[times_to_end_0_after_10 > 0])

mean1 <- mean(log(df1$time))
mean2 <- mean(log(df2$time))

p1 <- ggplot(df1, aes(x = log(time))) +
    geom_density() +
    geom_vline(aes(xintercept = exp(mean1)), color = "red", linetype = "dotted") +
    labs(title = paste("Density of Times to Reach 0 - Mean Time:", round(exp(mean1), 2)), x = "Time", y = "Density") + theme_bw()

p2 <- ggplot(df2, aes(x = log(time))) +
    geom_density() +
    geom_vline(aes(xintercept = exp(mean2)), color = "red", linetype = "dotted") +
    labs(title = paste("Density of Times to Reach 0 After Reaching 10 - Mean Time:", round(exp(mean2), 2)), x = "Time", y = "Density") + theme_bw() 

plot_df_0 <- do.call(rbind, paths_0)

p3 <- ggplot(plot_df_0, aes(x = log(time), y = y, group = sim)) +
    geom_line() +
    labs(title = "Paths of First Simulation", x = "Time", y = "Y") +
    theme_bw()

plot_df_0_after_10 <- do.call(rbind, paths_0_after_10)

p4 <- ggplot(plot_df_0_after_10, aes(x = log(time), y = y, group = sim)) +
    geom_line() +
    labs(title = "Paths of Second Simulation", x = "Time", y = "Y") +
    theme_bw()

grid.arrange(p1, p2, p3, p4, ncol = 2)

效率优化方法 #

1. 批量生成随机步长，减少采样开销 #

每次调用sample都会产生额外性能损耗，预先生成大量随机步长数组，模拟时直接取用，用完再补充：

# 预先生成1e6个步长
large_steps <- sample(c(-1, 1), 1e6, replace = TRUE, prob = c(0.5, 0.5))
step_ptr <- 1

# 模拟循环内直接调用预先生成的步长
while(y > 0) {
    y <- y + large_steps[step_ptr]
    step_ptr <- step_ptr + 1
    time <- time + 1
    # 步长耗尽时重新生成
    if(step_ptr > length(large_steps)) {
        large_steps <- sample(c(-1, 1), 1e6, replace = TRUE)
        step_ptr <- 1
    }
    if(y == 0) break
}

2. 避免动态扩展向量 #

c(path_0, y)会频繁重新分配内存，预先给路径向量分配足够大的初始空间，模拟结束后截断到实际长度：

max_steps <- 1e5  # 预设足够大的步数上限
path_0 <- numeric(max_steps)
path_0[1] <- 5
y <- 5
time <- 0
idx <- 2

while(y > 0 && idx <= max_steps) {
    y <- y + large_steps[step_ptr]
    step_ptr <- step_ptr + 1
    time <- time + 1
    path_0[idx] <- y
    idx <- idx + 1
    if(y == 0) break
}
# 截断到实际路径长度
path_0 <- path_0[1:(idx-1)]

3. 移除冗余操作，简化循环逻辑 #

删除循环内的print进度语句（仅在必要时批量输出），减少循环内的条件判断嵌套，提升单循环执行效率。

4. 采用解析解替代模拟 #

对于对称一维随机游走，这两个场景存在精确解析结论，无需模拟：

• 场景1：从位置5出发首次到达0的期望时间为无穷大，因为随机游走会无限徘徊，样本中会出现极长时间的路径，导致模拟无法收敛。
• 场景2：从位置5出发，先到达10再到达0的期望时间同样为无穷大，原因同上。

5. 并行化模拟 #

若坚持模拟，可使用R的并行包将模拟任务分配到多个核心运行：

library(parallel)
n_cores <- detectCores() - 1
cl <- makeCluster(n_cores)

# 定义场景1的单模拟函数
sim_scenario1 <- function() {
    y <- 5
    time <- 0
    while(y > 0) {
        y <- y + sample(c(-1,1), 1)
        time <- time + 1
        if(y == 0) break
    }
    time
}

# 并行运行100次模拟
times_to_end_0 <- parSapply(cl, 1:100, function(x) sim_scenario1())

stopCluster(cl)

内容的提问来源于stack exchange，提问作者stats_noob