R语言拒绝采样直方图无法绘制，请求代码修复

问题根源 #

1. 目标PDF逻辑矛盾：你描述的目标PDF仅在0≤x≤π/2时有值，但代码里的f函数额外定义了π/2 < x ≤ π的分支；更关键的是——提议PDFg(x)仅在π/2≤x≤π有值，和你描述的目标PDF支撑集几乎无交集，这种情况下拒绝采样根本无法获取目标分布的样本。我默认你是笔误，目标PDF应为0到π的三角形分布（和你代码原本的f逻辑一致），否则采样逻辑完全不成立。
2. M值计算错误：直接调用f(x)/g(x)时x未定义，且f、g是标量函数，无法处理向量输入；同时g(x)在0≤x<π/2时为0，会触发除以0的NaN报错。
3. 采样函数逻辑崩溃：接受样本时未递增计数变量i，会导致无限循环；return(samples)放在while循环内部，第一次循环就直接返回，无法生成指定数量的样本；当x落在0≤x<π/2时，g(x)=0会让接受概率计算出Inf，引发逻辑判断错误。

修复后的完整代码 #

# 目标PDF：0≤x≤π的三角形分布，峰值在π/2
f <- function(x) {
  ifelse(x >= 0 & x <= pi/2, (4/pi^2)*x,
         ifelse(x > pi/2 & x <= pi, (4/pi^2)*(pi - x), 0))
}

# 提议PDF：π/2≤x≤π时为sin(x)/2，其余为0
g <- function(x) {
  ifelse(x >= pi/2 & x <= pi, sin(x)/2, 0)
}

# 计算M：寻找f(x)/g(x)在重叠区间(π/2, π)的最大值
library(stats)
ratio_fun <- function(x) {
  if (x <= pi/2 | x >= pi) return(-Inf)
  f(x)/g(x)
}
# 用优化函数找最大值，乘1.01避免数值误差导致的采样失败
M <- optimize(ratio_fun, interval = c(pi/2, pi), maximum = TRUE)$objective * 1.01

# 修复后的拒绝采样函数
rejection_sampling <- function(n) {
  samples <- numeric(n)
  i <- 1
  while (i <= n) {
    # 直接从提议分布g(x)采样（用逆CDF方法，比均匀采样效率更高）
    u_g <- runif(1)
    x <- acos(1 - 2*u_g)
    
    # 防止数值误差导致x超出有效范围
    if (x < pi/2 | x > pi) next
    
    u <- runif(1)
    accept_prob <- f(x)/(M*g(x))
    # 检查接受概率有效性，符合条件则保留样本并递增计数
    if (!is.na(accept_prob) && u <= accept_prob) {
      samples[i] <- x
      i <- i + 1
    }
  }
  return(samples)
}

# 生成样本并绘制直方图
samples <- rejection_sampling(1000)
hist(samples, breaks=50, main="拒绝采样样本直方图", xlab="x")

核心修复说明 #

• 对齐分布支撑集：将目标PDF修正为和你代码逻辑一致的三角形分布，确保与提议分布有重叠区间，这是拒绝采样能工作的前提。
• 向量化函数改造：用ifelse替代if，让函数支持向量输入，方便后续最大值计算。
• 正确计算M值：使用optimize函数在有效区间内寻找比值最大值，避免除以0和未定义变量的问题。
• 提升采样效率：直接通过提议分布的逆CDF采样，无需从0到π盲目尝试，减少无效计算。
• 修复循环逻辑：仅在接受样本时递增计数变量，将return移至循环外部，确保生成足够数量的样本；增加NaN检查，避免逻辑判断报错。

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