如何修改ggplot与zoo代码实现多组OU过程同图对比（匹配base R效果）

问题背景 #

制作随机方法演示文稿时，需要在同一张图中对比多组参数不同的OU过程。已通过base R实现该效果，但现有ggplot代码会生成三张独立图表，需修改代码实现单图多组对比的可视化效果。

Base R实现代码 #

# Calculate the feature values according to an OU process
sim.ou <- function(theta, alpha, sigma, start_value) {
  x <- numeric(100)     # Time steps
  x[1] <- start_value   # Set the initial value
  for (i in 1:99) {
                        # OU simulation
    x[i + 1] <- x[i] - alpha * (x[i] - theta) + rnorm(n = 1, mean = 0, sd = sigma)
  }
  return(x)
}

# Plot the OU process
plotOU_diffusion <- function(sigma, theta, alpha, start_value, ...) {
  Xou <- replicate(1000, sim.ou(theta, alpha, sigma, start_value))
  matplot(Xou, type = "l", ylab = "", xlab = "", main = "", ...)
}

plotOU_diffusion(sigma=1,theta=1,alpha=1,start_value=1,col="#009E73",ylim=c(-7,7))
par(new=TRUE) 
plotOU_diffusion(sigma=0.5,theta=0.5,alpha=0.5,start_value=7,col="#56B4E9",ylim=c(-7,7))
par(new=TRUE)
plotOU_diffusion(sigma=0.1,theta=0.1,alpha=0.1,start_value=-4,col="#E69F00",ylim=c(-7,7))

Base R生成的效果 #

现有ggplot+zoo代码（生成三张独立图） #

library(tidyverse)
library(zoo)

# Calculate the feature values according to an OU process
sim.ou <- function(theta, alpha, sigma, start_value) {
  x <- numeric(150)     # Time steps
  x[1] <- start_value   # Set the initial value
  for (i in 1:149) {
                        # OU simulation
    x[i + 1] <- x[i] - alpha * (x[i] - theta) + rnorm(n = 1, mean = 0, sd = sigma)
  }
  return(x)
}

# Plot the OU process
plotOU_diffusion <- function(sigma, theta, alpha, start_value, num_walks, col, tit = NULL, xlabs = NULL, ylabs = NULL) {
  # Create a matrix of OU values
  Xou <- replicate(num_walks, sim.ou(theta, alpha, sigma, start_value))
  # Plot Xou object
  autoplot(zoo(Xou), facet = NULL, alpha = 0.5) + 
    theme_classic() + 
    geom_hline(yintercept = theta, linetype = "longdash", color = "black") +
    theme(legend.position = "none") +
    labs(title = tit, x = xlabs, y = ylabs) +
    scale_color_manual(values = rep(col,num_walks)) +
    scale_y_continuous(limits = c(-5, 5), expand = c(0,0)) +
    scale_x_continuous(limits = c(0, 150), expand = c(0,0))
}

plotOU_diffusion(sigma=1,theta=1,alpha=1,start_value=1,num_walks=50,col="#009E73")
plotOU_diffusion(sigma=0.5,theta=0.5,alpha=0.5,start_value=7,num_walks=50,col="#56B4E9")
plotOU_diffusion(sigma=0.1,theta=0.1,alpha=0.1,start_value=-4,num_walks=50,col="#E69F00")

现有ggplot生成的独立图表 #

修改后的ggplot+zoo代码（单图多组对比） #

核心思路是先将所有组的模拟数据整理为tidy格式的统一数据框，再通过一次ggplot调用完成所有曲线的绘制：

library(tidyverse)
library(zoo)

# 定义OU过程模拟函数（优化为可指定时间步数）
sim.ou <- function(theta, alpha, sigma, start_value, n_steps = 150) {
  x <- numeric(n_steps)
  x[1] <- start_value
  for (i in 1:(n_steps - 1)) {
    x[i + 1] <- x[i] - alpha * (x[i] - theta) + rnorm(n = 1, mean = 0, sd = sigma)
  }
  return(x)
}

# 整理所有组的参数信息
ou_groups <- tibble(
  group_id = factor(c("Group 1", "Group 2", "Group 3")),
  sigma = c(1, 0.5, 0.1),
  theta = c(1, 0.5, 0.1),
  alpha = c(1, 0.5, 0.1),
  start_val = c(1, 7, -4),
  line_color = c("#009E73", "#56B4E9", "#E69F00"),
  num_walks = c(50, 50, 50)
)

# 批量生成并整理所有模拟数据为长格式
ou_tidy_data <- ou_groups %>%
  rowwise() %>%
  mutate(
    # 生成当前组的所有路径
    paths = list(replicate(num_walks, sim.ou(theta, alpha, sigma, start_val))) %>%
      # 转换为数据框并添加时间步、路径ID列
      map(~ as_tibble(.x) %>%
            mutate(time_step = row_number()) %>%
            pivot_longer(-time_step, names_to = "path_id", values_to = "value"))
  ) %>%
  unnest(paths)

# 绘制单图对比所有组的OU过程
ggplot(ou_tidy_data, aes(x = time_step, y = value, color = group_id)) +
  # 绘制所有路径曲线，设置透明度避免重叠过密
  geom_line(alpha = 0.5) +
  # 添加每组的theta水平参考线
  geom_hline(aes(yintercept = theta), data = ou_groups, linetype = "longdash", color = "black") +
  # 使用经典主题
  theme_classic() +
  # 指定各组的颜色
  scale_color_manual(values = ou_groups$line_color) +
  # 设置坐标轴范围（对齐base R的效果）
  scale_y_continuous(limits = c(-7, 7), expand = c(0, 0)) +
  scale_x_continuous(limits = c(1, 150), expand = c(0, 0)) +
  # 隐藏坐标轴标签（对齐base R效果）
  labs(x = "", y = "") +
  # 隐藏图例（如需保留可删除此行）
  theme(legend.position = "none")

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