如何在R语言中从自定义PDF分布生成10000个随机样本？

针对特定分布：( f(x) = 3x^2 I(0 < x < 1) ) #

这个分布的累积分布函数（CDF）可通过积分推导：

( F(x) = \int_0^x 3t^2 dt = x^3 )，其中 ( 0 < x < 1 )

其逆CDF为 ( F^{-1}(u) = u^{1/3} )（( u ) 服从均匀分布 ( U(0,1) )）。利用逆变换法，生成10000个样本的代码如下：

# 设置随机种子保证结果可复现
set.seed(123)
# 生成10000个U(0,1)均匀样本
u_samples <- runif(10000)
# 通过逆CDF转换得到目标分布样本
target_samples <- u_samples^(1/3)

验证样本有效性 #

可以通过直方图对比理论概率密度函数（PDF）来验证样本是否符合分布：

hist(target_samples, breaks = 30, freq = FALSE, 
     main = "样本直方图与理论PDF", xlab = "x")
curve(3*x^2, from = 0, to = 1, col = "red", lwd = 2, add = TRUE)
legend("topright", legend = "理论PDF", col = "red", lwd = 2)

从自定义PDF生成随机样本的通用方法 #

1. 逆变换法（首选场景：逆CDF可求解） #

核心逻辑：利用均匀分布样本，通过目标分布的逆CDF转换得到符合要求的样本。

• 步骤：
  • 计算目标分布的CDF ( F(x) = \int_{-\infty}^x f(t)dt )
  • 求解逆函数 ( F^{-1}(u) )（( u \sim U(0,1) )）
  • 生成U(0,1)样本，代入逆函数得到目标样本
• 适用场景：CDF逆函数容易推导的情况（如本次示例、指数分布、均匀分布等）

2. 接受-拒绝法（适用场景：逆变换不可用） #

核心逻辑：通过一个易采样的提议分布，筛选出符合目标分布的样本。

• 步骤：
  • 选择提议分布 ( g(x) )（如均匀、正态分布），找到常数 ( M ) 使得 ( f(x) \leq M \cdot g(x) ) 对所有x成立
  • 生成提议样本 ( y \sim g(x) ) 和均匀样本 ( u \sim U(0,1) )
  • 若 ( u \leq \frac{f(y)}{M \cdot g(y)} )，则接受y作为样本；否则拒绝并重复操作
• 示例代码（以目标PDF ( f(x)=4x^3 I(0<x<1) ) 为例）：

set.seed(123)
n <- 10000
samples <- numeric(n)
i <- 1
# 提议分布用U(0,1)，M取max(f(x)/g(x))=4
while(i <= n) {
  y <- runif(1)
  u <- runif(1)
  if(u <= (4*y^3)/(4*1)) {
    samples[i] <- y
    i <- i + 1
  }
}

3. Metropolis算法（适用场景：高维或PDF复杂，MCMC方法） #

核心逻辑：通过马尔可夫链逐步生成符合目标分布的样本，适合高维或无法直接计算归一化常数的情况。

• 步骤：
  • 初始化样本 ( x_0 )
  • 从提议分布（如正态分布 ( N(x_{t-1}, \sigma^2) )）生成候选样本 ( y )
  • 计算接受概率 ( \alpha = \min\left(1, \frac{f(y)}{f(x_{t-1})}\right) )
  • 生成 ( u \sim U(0,1) )，若 ( u \leq \alpha ) 则更新样本为y，否则保留原样本
  • 重复步骤2-4，丢弃前一段"burn-in"样本（消除初始值影响）
• 示例代码（针对本次目标分布）：

set.seed(123)
n <- 10000
burn_in <- 1000 # 丢弃前1000个burn-in样本
samples <- numeric(n + burn_in)
samples[1] <- 0.5 # 初始值
sigma <- 0.2 # 提议分布的标准差

for(t in 2:(n + burn_in)) {
  y <- rnorm(1, mean = samples[t-1], sd = sigma)
  # 确保候选样本在(0,1)范围内，否则目标PDF值为0，不接受
  if(y <= 0 || y >= 1) {
    samples[t] <- samples[t-1]
    next
  }
  # 计算接受概率（目标PDF的比值）
  alpha <- min(1, (3*y^2)/(3*samples[t-1]^2))
  u <- runif(1)
  samples[t] <- ifelse(u <= alpha, y, samples[t-1])
}
# 最终样本（丢弃burn-in）
final_samples <- samples[(burn_in + 1):(n + burn_in)]

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