在R中计算多项分布真实分布与MLE/贝叶斯估计分布的总变差距离

嘿，我来帮你搞定这个多项分布下MLE和贝叶斯估计量的对比任务，包括样本生成、估计量计算和总变差距离的统计，具体步骤和R代码如下：

你要完成的核心任务是：

• 从真实多项分布生成400个独立样本，每个样本对应30次试验的10类别计数结果
• 对每个样本计算10个概率参数的MLE和贝叶斯估计量
• 分别统计两种估计量与真实分布之间的总变差距离（TVD）

注意：rmultinom的prob参数要求和为1，你给出的权重c(5,7,10,8,14,10,15,12,10,9)总和刚好是100，直接归一化就能得到真实概率向量。

用R的rmultinom生成样本，我们将结果转置为每行对应一个样本的格式，方便后续处理：

# 定义真实概率（归一化权重）
true_prob <- c(5,7,10,8,14,10,15,12,10,9)/100
# 设置随机种子保证结果可复现
set.seed(123)
# 生成400个样本，每个样本是30次试验的类别计数
samples <- rmultinom(n = 400, size = 30, prob = true_prob)
# 转置矩阵：400行（样本）×10列（类别计数）
samples <- t(samples)

3.1 MLE估计量 #

多项分布的MLE非常直观：每个类别的计数除以总试验次数（30）即可得到概率估计：

# 批量计算所有样本的MLE：每行对应一个样本的10个概率估计
mle_estimates <- samples / 30

3.2 贝叶斯估计量（共轭Dirichlet先验） #

通常用Dirichlet共轭先验来计算贝叶斯估计量，最常用的是均匀先验Dirichlet(1,1,...,1)（对应拉普拉斯平滑），此时估计量公式为：
$$\hat{p}_i = \frac{n_i + 1}{30 + 10}$$
其中$n_i$是第i类的计数，分母是总试验次数加上10个先验参数的总和。

你也可以根据需求调整先验参数（比如Dirichlet(0.5,...,0.5)的Jeffreys先验），这里以均匀先验为例：

# 批量计算所有样本的贝叶斯估计量
bayes_estimates <- (samples + 1) / (30 + 10)

总变差距离的定义是两个概率向量$P$和$Q$的L1距离的一半：
$$TVD(P, Q) = \frac{1}{2} \sum_{i=1}^{k} |P_i - Q_i|$$
其中k是类别数（这里k=10）。

我们可以写一个函数批量计算每个样本的TVD：

# 定义TVD计算函数
compute_tvd <- function(est_prob, true_prob) {
  0.5 * sum(abs(est_prob - true_prob))
}

# 计算每个样本的MLE与真实分布的TVD
mle_tvd <- apply(mle_estimates, 1, compute_tvd, true_prob = true_prob)
# 计算每个样本的贝叶斯估计量与真实分布的TVD
bayes_tvd <- apply(bayes_estimates, 1, compute_tvd, true_prob = true_prob)

你可以对TVD结果做统计分析，或者绘制箱线图直观对比两种估计量的表现：

# 查看TVD的描述性统计
summary(mle_tvd)
summary(bayes_tvd)

# 绘制箱线图对比两种估计量的TVD分布
boxplot(mle_tvd, bayes_tvd, 
        names = c("MLE TVD", "Bayes TVD"),
        ylab = "Total Variation Distance",
        main = "TVD Comparison: MLE vs Bayesian Estimator")

通过这些结果就能清晰看到哪种估计量更接近真实分布。

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