如何在R中基于模拟结果计算P(log(Y) > sin(X))的概率近似值

实现原理 #

你提供的代码已经通过拒绝采样得到了服从目标联合分布的有效样本集，蒙特卡洛法估计概率的核心逻辑就是：满足条件log(Y) > sin(X)的样本数占总有效样本数的比例，即为所求概率的近似值。

具体代码实现 #

完整流程代码（包含模拟+概率计算） #

# 1. 生成有效样本（若你已经完成模拟可直接跳过该段）
# 可根据精度需求调大初始样本量n，样本量越大结果越准确
n <- 1e6 
sim_raw <- data.frame(
  x = runif(n, -1, 1),
  y = runif(n, 0, 1),
  z = runif(n, 0, 3/2)
)
# 过滤得到符合目标分布的有效样本
valid_idx <- with(sim_raw, 0 < y & x^2 + y^2 < 1 & z <= (3/2)*y) 
valid_sample <- sim_raw[valid_idx, ]

# 2. 计算概率近似值
# 逻辑向量的均值即为满足条件的样本占比，等价于sum(条件)/总样本数
prob_est <- mean(log(valid_sample$y) > sin(valid_sample$x))

# 输出结果
cat("P(log(Y) > sin(X))的近似值：", round(prob_est, 4), "\n")
cat("本次计算使用的有效样本量：", nrow(valid_sample), "\n")

可选：计算估计值的95%置信区间 #

如果需要评估估计结果的可靠性，可以添加以下代码计算置信区间：

n_valid <- nrow(valid_sample)
se <- sqrt(prob_est * (1 - prob_est) / n_valid)
ci_95 <- c(prob_est - 1.96 * se, prob_est + 1.96 * se)
cat("95%置信区间：(", round(ci_95[1],4), ",", round(ci_95[2],4), ")\n")

注意事项 #

• 初始样本量n越大，有效样本量越多，估计结果的精度越高。如果n取1e4的话有效样本大概在2500左右，精度较低，建议调至1e5或1e6使用。
• 代码中直接对逻辑向量取均值是R的简化写法，效果和sum(log(valid_sample$y) > sin(valid_sample$x))/nrow(valid_sample)完全一致。

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