基于混合分布从密度函数生成伪随机数的技术求助

问题需求 #

生成4900个与以下密度函数成比例的伪随机数，存储至变量X3：

• 当x < -6时，f(x) = 0
• 当-6 ≤ x < 0时，f(x) = 4
• 当x ≥ 0时，f(x) = 3*exp(-2.1x) + 2*exp(-1.4x)

课堂参考案例 #

参考任务 #

基于混合分布生成10000个与如下密度函数成比例的数值，近似计算随机变量均值：

• 当x ≥ -1时，f(x) = exp(-2x²) + 0.5*exp(-2x)
• 当x < -1时，f(x) = 0

参考R代码及注释 #

# 定义成比例的目标密度函数
f <- function(x){
  return((exp(-2*x^2)+0.5*exp(-2*x))*(x>= -1))
}
# 计算正态分布N(0,0.5²)在x≥-1区间的累积概率上下限
F1_a <- pnorm(-1,sd=0.5) 
F1_b <- 1
# 绘制目标密度曲线
curve(f,-2,10)

# 定义第一个混合分量的条件密度：N(0,0.5²)在x≥-1的条件分布
d1 <- function(x){
  return(dnorm(x,sd=0.5)/(F1_b-F1_a)*(x>=-1))
}
# 定义第二个混合分量的密度：平移后的指数分布
d2 <- function(x){return(dexp(x+1,rate=2))}
# 计算目标密度的总积分（用于归一化）
int_f <- (sqrt(2*pi)*0.5)*(F1_b-F1_a)+exp(2)/4 
# 计算两个混合分量的权重
p1 <- (sqrt(2*pi)*0.5)*(F1_b-F1_a)/int_f
p2 <- exp(2)/4/int_f
# 验证混合分解正确性（曲线应趋近于0）
curve(f(x)/int_f-p1*d1(x)-p2*d2(x),-2,5)
# 注：微小差异来自浮点计算误差

# 生成10000个样本
n<-10000
# 按权重选择样本所属分量
milline <- sample.int(2,n,replace=TRUE,prob = c(p1,p2))
tulem <- rep(NA,n)
mitu <- table(milline)
# 从第一个分量抽样：正态分布条件抽样
tulem[milline==1]<- qnorm(runif(mitu[1],F1_a,F1_b),sd=0.5)
# 从第二个分量抽样：指数分布平移变换
tulem[milline==2]<- rexp(mitu[2],2)-1
# 绘制直方图与真实密度对比
hist(tulem,30,probability =TRUE)
curve(f(x)/int_f,-1,5,add=TRUE)
# 计算样本均值
mean(tulem)

当前问题的实现方案 #

核心逻辑 #

目标密度可拆分为两个有效混合分量（x < -6概率为0，无需处理）：

1. 分量2：-6 ≤ x < 0的均匀分布（密度为常数4）
2. 分量3：x ≥ 0的两个指数分布混合（3*exp(-2.1x)和2*exp(-1.4x)）

实现步骤：

1. 计算各分量的积分占总积分的比例（即权重）
2. 按权重随机分配每个样本到对应分量
3. 针对分量特性生成伪随机数

R代码实现 #

# 定义成比例的目标密度函数
target_f <- function(x) {
  ifelse(x < -6, 0,
         ifelse(x >= -6 & x < 0, 4,
                3*exp(-2.1*x) + 2*exp(-1.4*x)))
}

# 计算各分量的积分
# 分量2的积分：[-6,0]区间上常数4的积分
int_comp2 <- 4 * (0 - (-6))  # 结果为24
# 分量3的积分：两个指数项的积分和（∫0到∞ exp(-λx)dx = 1/λ）
int_comp3 <- 3*(1/2.1) + 2*(1/1.4)  # 约2.8572
# 总积分
int_total <- int_comp2 + int_comp3

# 计算分量权重
p2 <- int_comp2 / int_total
p3 <- int_comp3 / int_total

# 生成4900个样本
n <- 4900
X3 <- rep(NA, n)
# 随机选择每个样本所属分量
component <- sample.int(2, n, replace = TRUE, prob = c(p2, p3))
count_comp2 <- sum(component == 1)
count_comp3 <- sum(component == 2)

# 分量2：生成[-6,0]上的均匀分布随机数
X3[component == 1] <- runif(count_comp2, min = -6, max = 0)

# 分量3：处理两个子指数分布的混合
# 计算子分量权重
w3a <- (3/2.1) / int_comp3  # 对应3*exp(-2.1x)的权重
w3b <- (2/1.4) / int_comp3  # 对应2*exp(-1.4x)的权重
# 为分量3的样本选择子分量
sub_component <- sample.int(2, count_comp3, replace = TRUE, prob = c(w3a, w3b))
count_sub3a <- sum(sub_component == 1)
count_sub3b <- sum(sub_component == 2)
# 生成子分量样本：指数分布随机数
X3[component == 2 & sub_component == 1] <- rexp(count_sub3a, rate = 2.1)
X3[component == 2 & sub_component == 2] <- rexp(count_sub3b, rate = 1.4)

# 验证：绘制样本直方图与归一化后的目标密度对比
hist(X3, breaks = 50, probability = TRUE, main = "样本分布 vs 目标密度")
curve(target_f(x)/int_total, from = -6, to = 5, add = TRUE, col = "red", lwd = 2)

关键说明 #

• 分量2的均匀分布直接用runif生成区间内随机数，符合常数密度的特性
• 分量3的两个指数项，因密度形式与指数分布Exp(λ)（密度λexp(-λx)）成比例，直接用rexp(rate=λ)生成即可，再按各自系数的积分占比分配权重
• 最终变量X3即为符合要求的4900个伪随机数

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