为何下述R代码生成的变量x是五个伽马分布的混合样本？

这段代码是伽马混合分布的标准向量化采样实现，你提到的混合系数$\theta_k$已经在代码里明确定义，没有缺失，逻辑对应关系如下：

有限混合分布的采样逻辑 #

有限混合分布的通用采样流程本来就分为两步：

• 第一步：对每个待生成的样本，按照混合系数$\theta_k$的概率分布，抽样确定该样本归属的混合分量编号k
• 第二步：根据每个样本所属分量的对应分布参数，从分量分布中生成该样本的取值
  所有样本生成后拼接得到的整体样本，就是对应混合模型的采样结果。

代码和模型的对应关系 #

逐行对应你的代码看：

1. 混合系数设定：sample(1:5, size=n, replace=TRUE, prob=(1:5)/15)里的prob=(1:5)/15就是你要的混合系数$\theta_k$，计算可得5个分量的权重分别为$\theta_1=1/15、\theta_2=2/15、\theta_3=3/15、\theta_4=4/15、\theta_5=5/15$，权重和为1，完全符合混合系数的要求。这行代码执行后，长度为5000的向量k里，每个位置的取值就是对应样本归属的分量编号，每个编号被抽中的概率严格等于对应混合系数。
2. 分量参数设定：rate <- 1/k是向量化的参数赋值，对应第k个伽马分量的率参数为1/k，所有分量的形状参数统一为3，即5个混合分量分别为$Gamma(shape=3, rate=1/1)$、$Gamma(shape=3, rate=1/2)$……$Gamma(shape=3, rate=1/5)$，和模型定义的分量完全匹配。
3. 分量采样：R的rgamma函数支持向量化传参，执行x <- rgamma(n, shape=3, rate=rate)时，会对第i个样本使用rate[i]作为率参数生成伽马随机数，等价于每个样本从自己被分配到的分量中抽样取值。

这种写法和“先按混合系数计算每个分量要抽取的样本量，再循环对每个分量抽样后拼接”的传统循环写法结果完全等价，只是利用R的向量化特性省略了循环步骤，运行效率更高。

你可以自行做简单验证：抽样完成后统计k向量中1-5的取值占比，结果会非常接近你设定的$\theta_k$；按k的取值分组统计x的均值、方差，也会和对应伽马分量的理论矩匹配。

n <- 5000
k <- sample(1:5, size=n, replace=TRUE, prob=(1:5)/15)
rate <- 1/k
x <- rgamma(n, shape=3, rate=rate)
# 验证分量占比
prop.table(table(k))
# 验证分组矩
tapply(x, k, mean)

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