r - 为什么 `x` 现在是以下 R 代码中这五个 Gamma 分布的混合？

Question

我尝试从以下 Gamma 分布的混合模型中进行采样:

enter image description here

R代码如下:

该算法可以转化为矢量化方法。

第 1 步:生成随机样本 k_1,...,k_n向量中的整数 k ，其中P(k)=θ_k ，k=1，...，5。

第2步:设置rate=1/k。

n <- 5000
k <- sample(1:5, size=n, replace=TRUE, prob=(1:5)/15)
rate <- 1/k
x <- rgamma(n, shape=3, rate=rate)

我的问题是为什么x现在是这五个 Gamma 分布的混合吗？在混合模型的表达中，似乎还需要系数theta_k ？

Answer 1

这里有两种方法可以将 Gamma 混合分布的样本与预期的混合密度进行比较。这应该有助于理解 F_X 如何是 Gamma 分布混合的累积分布函数。

# Fix random seed for reproducibility
set.seed(2022)

# Sample data
n <- 100000
X <- unlist(lapply(1:5, function(j) rgamma(n * j, shape = 3, rate = 1 / j)))

# Weighted Gamma mixture density
dmix <- function(x)  sapply(
    x,
    function(val) sum((1:5) / 15 * dgamma(val, shape = 3, rate = 1 / (1:5))))
library(tibble)
data = tibble(x = seq(0, ceiling(max(X)), length.out = 100), y = dmix(x))

# Plot histogram of samples and compare with density 
library(ggplot2)
ggplot(data.frame(x = X), aes(x)) +
    geom_histogram(aes(y = ..density..), bins = 200) + 
    geom_line(data = data, aes(x, y))

评论:

1. 采样时，我们根据权重调整每个 Gamma 分布的样本数量；在本例中，权重仅为 1:5。
2. 我们使用 aes(y = ..密度..) 将直方图表示为正确归一化的密度，以便我们可以将值与混合密度 dmix 进行比较。