在R中基于频率表从截断多元正态分布抽样的实现方案

需求说明 #

需要从指定参数的二元正态分布中抽取样本，参数如下：

mu = c(32, 170) # 第一个元素为年龄均值，第二个为身高均值
sigma = matrix(c(9, 10,
                10, 100), nrow = 2)

要求样本年龄符合以下频率分布：

原尝试的循环抽样代码存在筛选条件错误（误将样本数量作为区间阈值），且反复抽样筛选的方式效率极低，甚至可能出现无限循环。

高效解决方案：利用条件正态分布直接抽样 #

不需要反复抽样筛选，而是基于条件多元正态分布的性质，直接生成符合要求的样本，步骤如下：

1. 推导条件分布参数 #

对于二元正态分布N(μ, Σ)：

• 年龄（X₁）的边缘分布：N(μ₁=32, σ₁₁=9)
• 当已知X₁=x₁时，身高（X₂）的条件分布为：
  • 条件均值：$\mu_{2|1} = \mu_2 + \frac{\sigma_{12}}{\sigma_{11}}(x_1 - \mu_1)$
  • 条件方差：$\sigma_{2|1}^2 = \sigma_{22} - \frac{\sigma_{12}^2}{\sigma_{11}}$

代入参数计算：

mu1 <- mu[1]
mu2 <- mu[2]
sigma11 <- sigma[1,1]
sigma12 <- sigma[1,2]
sigma22 <- sigma[2,2]

# 条件分布参数
cond_mean <- function(x1) mu2 + (sigma12 / sigma11) * (x1 - mu1)
cond_var <- sigma22 - (sigma12^2) / sigma11

2. 生成符合区间要求的年龄样本 #

使用截断正态分布生成每个区间内的年龄样本，这里用truncnorm包的rtruncnorm函数实现：

# 安装并加载truncnorm包（如果未安装）
# install.packages("truncnorm")
library(truncnorm)

freq_table <- data.frame(age_low = c(25, 30, 35), n = c(20, 43, 66))

# 生成年龄样本
age_samples <- pmap(freq_table, function(age_low, n) {
  rtruncnorm(n = n, a = age_low, b = age_low + 1, mean = mu1, sd = sqrt(sigma11))
}) %>% unlist()

3. 生成对应的身高样本 #

基于每个年龄样本值，生成对应的身高样本：

height_samples <- rnorm(length(age_samples), mean = cond_mean(age_samples), sd = sqrt(cond_var))

4. 组合成最终样本矩阵 #

final_samples <- cbind(age_samples, height_samples)
head(final_samples)

用dplyr实现的简洁版本 #

如果习惯用管道式写法，可以用dplyr+purrr组合实现：

library(dplyr)
library(purrr)
library(truncnorm)

freq_table <- tibble(age_low = c(25, 30, 35), n = c(20, 43, 66))

final_samples <- freq_table %>%
  mutate(
    age = map2(age_low, n, ~rtruncnorm(.y, a = .x, b = .x + 1, mean = mu1, sd = sqrt(sigma11))),
    height = map(age, ~rnorm(length(.x), mean = cond_mean(.x), sd = sqrt(cond_var)))
  ) %>%
  unnest(c(age, height)) %>%
  as.matrix()

原代码问题说明 #

原代码的核心错误是筛选条件写错：

# 原错误条件：把样本数量n当成了年龄区间
sample[sample[,1] >= frqtable[row, 'n'] & sample[,1] < frqtable[row, 'n'] + 1,]
# 正确应该是用age列：
sample[sample[,1] >= frqtable[row, 'age'] & sample[,1] < frqtable[row, 'age'] + 1,]

即使修正条件，反复抽样筛选的方式效率极低，尤其是当目标区间的概率较小时，会浪费大量内存和时间，甚至陷入无限循环。而基于条件分布的方法直接生成符合要求的样本，完全避免了筛选步骤，效率提升显著。

内容的提问来源于stack exchange，提问作者Aku-Ville Lehtimäki