我的答案是关于如何使用基数 R 来计算。
我们有两个多项式参数向量，θ 和 η。总变化距离等价于P_θ(E) - P_η(E)，其中E={ω | P_θ({ω})>P_η({ω})}，ω是样本数的向量。
我知道在基础 R 中评估 P(E) 的两种方法。一种是非常简单的基于模拟的方法。另一个根据计数的线性组合（近似正态分布）来重构问题，并使用pnorm 函数。
基于模拟的方法
您模拟来自每个分布的样本，使用概率质量函数检查它们是否在 E 中，并计算它们的频率。我将在这里举一个例子。我们将根据您的问题假设真实分布：
unnormalized.true <- c(5,7,10,8,14,10,15,12,10,9)
true <- unnormalized.true / sum(unnormalized.true)
我们将抽取样本并使用贝叶斯估计器估计新分布：
set.seed(921)
result <- as.vector(rmultinom(1, size = 30, prob = true))
result
## [1] 3 6 2 0 5 3 3 4 1 3
dirichlet <- (result+1)/(30+length(true))
计算真实分布下E的概率：
set.seed(939)
true.dist <- rmultinom(10^6, 30, true)
p.true.e <- mean(apply(true.dist, 2, function(x)
dmultinom(x, 30, true) - dmultinom(x, 30, dirichlet) > 0))
在贝叶斯估计器的估计分布下计算E的概率：
dirichlet.dist <- rmultinom(10^6, 30, dirichlet)
p.dirichlet.e <- mean(apply(dirichlet.dist, 2, function(x)
dmultinom(x, 30, true) - dmultinom(x, 30, dirichlet) > 0))
我们可以减去得到总变化距离。
p.true.e - p.dirichlet.e
## [1] 0.83737
用最大似然估计重复这个，我们得到估计量的比较。
mle <- result/30
mle.dist <- rmultinom(10^6, 30, mle)
p.true.e2 <- mean(apply(true.dist, 2, function(x)
dmultinom(x, 30, true) - dmultinom(x, 30, mle) > 0))
p.mle.e2 <- mean(apply(mle.dist, 2, function(x)
dmultinom(x, 30, true) - dmultinom(x, 30, mle) > 0))
p.true.e2 - p.mle.e2
## [1] 0.968301
（已编辑以修复一个严重错误。之前我在与 MLE 的比较中重新使用了p.true.e。我忘记了事件 E 是根据估计分布来定义的。）
正态近似
我认为这种方法实际上比基于模拟的方法更准确，尽管是正态近似。正如您将看到的，我们没有对多项式计数进行正态近似，这对于 n=30 不太可能准确。我们正在对这些计数的线性组合进行正态近似，这接近于正态。这种方法的缺点是它不能处理估计分布中的零概率。这是一个真正的问题，因为对我来说，优雅地处理零点是使用总变化距离而不是 Kullback-Leibler 散度的一部分。但它就在这里。
以下推导产生E的重述：
定义
其中 N_i 是多项式样本的一个单元格，并且
那么，E就是L>0的事件。
我们遇到零概率问题的原因是它导致 λ_i 之一是无限的。
我想验证 L 接近正态分布，在前面的示例中。我将通过使用前面的多项式模拟从 L 的分布中获取样本来做到这一点：
lambda <- log(true/dirichlet)
L.true.dist <- apply(true.dist, 2, function(x) sum(lambda*x))
L.dirichlet.dist <- apply(dirichlet.dist, 2, function(x) sum(lambda*x))
请注意，我正在比较真实分布和贝叶斯估计分布。我不能用 MLE 做那个，因为我的样本计数为零。
绘制 L 的分布并与正态拟合进行比较：
par(mfrow=c(1,2))
L.true.dist.hist <- hist(L.true.dist)
L.true.dist.fit <- function(x)
length(L.true.dist) * diff(L.true.dist.hist$breaks)[1] *
dnorm(x, mean(L.true.dist), sd=sd(L.true.dist))
curve(L.true.dist.fit, add=TRUE, n=1000, col='red')
L.dirichlet.dist.hist <- hist(L.dirichlet.dist)
L.dirichlet.dist.fit <- function(x)
length(L.dirichlet.dist) * diff(L.dirichlet.dist.hist$breaks)[1] *
dnorm(x, mean(L.dirichlet.dist), sd=sd(L.dirichlet.dist))
curve(L.dirichlet.dist.fit, add=TRUE, n=1000, col='red')
par(mfrow=c(1,1))
L 的分布看起来很正常。因此，我们可以使用pnorm，而不是使用模拟。但是，我们需要计算 L 的均值和标准差。这可以按如下方式完成。
L 的均值是
其中p_i 是细胞i 在分布p 中的细胞概率。方差是
在哪里
是多项分布的协方差矩阵。我将为此示例计算这些矩，并根据模拟中的经验矩检查它们。一、对于真实分布下L的分布：
n <- 30
k <- length(true)
mean.L.true <- sum(lambda * n * true)
# Did we get the mean right?
c(mean.L.true, mean(L.true.dist))
## [1] 3.873509 3.875547
# Covariance matrix assuming the true distribution
sigma.true <- outer(1:k, 1:k, function(i,j)
ifelse(i==j, n*true[i]*(1-true[i]), -n*true[i]*true[j]))
var.L.true <- t(lambda) %*% sigma.true %*% lambda
# Did we get the standard deviation right?
c(sqrt(var.L.true), sd(L.true.dist))
## [1] 2.777787 2.776945
那么，L在分布的贝叶斯估计下的均值和方差：
mean.L.dirichlet <- sum(lambda * n * dirichlet)
# Did we get the mean right?
c(mean.L.dirichlet, mean(L.dirichlet.dist))
## [1] -3.893836 -3.895983
# Covariance matrix assuming the estimated distribution
sigma.dirichlet <- outer(1:k, 1:k, function(i,j)
ifelse(i==j, n*dirichlet[i]*(1-dirichlet[i]), -n*dirichlet[i]*dirichlet[j]))
var.L.dirichlet <- t(lambda) %*% sigma.dirichlet %*% lambda
# Did we get the standard deviation right?
c(sqrt(var.L.dirichlet), sd(L.dirichlet.dist))
## [1] 2.796348 2.793421
有了这些，我们就可以用pnorm计算总变异距离：
pnorm(0, mean.L.true, sd=sqrt(var.L.true), lower.tail=FALSE) -
pnorm(0, mean.L.dirichlet, sd=sqrt(var.L.true), lower.tail=FALSE)
## [1] 0.8379193
# Previous result was 0.83737
我们得到三位数的模拟结果。
不过，我不知道有什么简单的方法可以扩展正态近似方法来处理零概率。我有一个想法，但我在试图计算计数的协方差矩阵时遇到了困难，条件是计数为 0 的特定单元格。如果您认为可以有所作为，我可以分享我的进展。