基于概率密度函数的While循环代码故障：x_accept为空问题求助

问题背景 #

需求：编写While循环，循环执行至n等于105时停止；每次生成15-33区间的随机值x1，同时生成0-1均匀分布的随机变量num，当x1满足条件时将其存入x_accept。

原运行代码如下：

library(ggplot2)
dev.new()
n = 0
set.seed(1684)
x = seq(15, 33, by = 0.1)
f <- function(x) {
  out <- ifelse(
    x < 15 | 33 < x,
    0,
    ifelse(
      15 <= x & x <= 24,
      (2*(x-15))/((33-15)*(24-15)),
      ifelse(
        24 < x & x <= 33,
        (2*(33-x))/((33-15)*(33-24)),
        NA_real_
      )))
  if (any((is.na(out) | is.nan(out)) & (!is.na(x) & !is.nan(x)))) {
    warning("f(x) undefined for some input values")
  }
  out
}

while (n != 105) {
  n = n + 1
  x1 = runif(1, min = 15 , max = 33)
  num = runif(1, min = 0 , max = 1)
  if ((num < (f(x1)/2/(33-15))) && (num == (18*(f(x1)/2)))) {
    
    x_accept = c(x_accept, x1)
  }
}

运行后发现x_accept始终为空，以下是问题排查与修复方案：

问题原因分析 #

1. 逻辑条件完全无法触发
   原代码的if判断使用&&同时要求两个条件成立：num < ...和num == ...。但num是连续型均匀分布随机数，等于某个精确浮点值的概率为0，这直接导致条件永远无法满足，x_accept永远不会被赋值。

2. 拒绝采样的接受概率公式错误
   目标分布f(x)是一个三角形分布，在x=24处取得最大值：

   f(24) = (2*(24-15))/((33-15)*(24-15)) = 18/(18*9) = 1/9 ≈ 0.1111
   

   拒绝采样的核心逻辑是num < f(x1)/M（M为目标分布的最大值），原代码中的f(x1)/2/(33-15)计算错误，实际应为f(x1)/M，其中M=1/9。

3. 变量未初始化
   x_accept在首次赋值前未定义，即使条件触发也会抛出变量未找到的错误。

修复后的代码 #

library(ggplot2)
dev.new()
n = 0
set.seed(1684)
x = seq(15, 33, by = 0.1)
f <- function(x) {
  out <- ifelse(
    x < 15 | 33 < x,
    0,
    ifelse(
      15 <= x & x <= 24,
      (2*(x-15))/((33-15)*(24-15)),
      ifelse(
        24 < x & x <= 33,
        (2*(33-x))/((33-15)*(33-24)),
        NA_real_
      )))
  if (any((is.na(out) | is.nan(out)) & (!is.na(x) & !is.nan(x)))) {
    warning("f(x) undefined for some input values")
  }
  out
}

# 初始化存储变量
x_accept = c()
# 计算目标分布的最大值M
M = 1/9

while (n != 105) {
  n = n + 1
  x1 = runif(1, min = 15 , max = 33)
  num = runif(1, min = 0 , max = 1)
  # 修正接受条件：仅判断num是否小于f(x1)/M
  if (num < f(x1)/M) {
    x_accept = c(x_accept, x1)
  }
}

# 查看采样结果
print(x_accept)
# 绘制采样结果与目标分布对比
ggplot(data.frame(x=x_accept), aes(x=x)) +
  geom_histogram(aes(y=after_stat(density)), bins=15, alpha=0.5) +
  stat_function(fun=f, color="red", linewidth=1)

验证说明 #

修复后代码会正确执行拒绝采样：

• 初始化x_accept避免未定义错误
• 修正接受条件为标准拒绝采样逻辑，确保符合目标三角形分布的样本被接受
• 循环执行105次后，x_accept会存储符合条件的采样值，可通过打印或绘图验证分布是否与f(x)一致

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