找到 r 中多项式分布之间的总变异距离
find total variation distance between multinomial distributions in r
我正在将贝叶斯估计量与多项式分布中的 MLE 进行比较。我正在使用
从特定的多项式分布中使用 rmultinom 绘制随机样本
rmultinom(400, size = 30, prob = c(5,7,10,8,14,10,15,12,10,9))
对于 400 个样本中的每一个，我计算了十个概率参数的 MLE 和贝叶斯估计量。我现在想在每种情况下找到真实分布与估计量定义的分布之间的总变异距离。
由于 30 号和 10 号箱有超过 2 亿种可能的排列方式，我认为使用理论定义不是一个好主意。
包 distrEx has a function "TotalVarDist()", but it can only be used with distributions defined in the distr package, and multinomial is not one of them. There are directions for defining them (see here and here) 但选项是通过显式列出支持来定义离散分布（同样，我认为这不是一个好的选择，因为支持的大小超过 200百万）或使用与创建 distr 包相同的方法从头开始，这超出了我目前的能力。
关于如何使用提到的包或以完全不同的方式执行此操作的任何想法？
我的回答是关于如何使用基数 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 的特定单元格为条件的计数的协方差矩阵时遇到了困难。如果你认为你可以有所作为，我可以分享我的进步。
我正在将贝叶斯估计量与多项式分布中的 MLE 进行比较。我正在使用
从特定的多项式分布中使用 rmultinom 绘制随机样本rmultinom(400, size = 30, prob = c(5,7,10,8,14,10,15,12,10,9))
对于 400 个样本中的每一个，我计算了十个概率参数的 MLE 和贝叶斯估计量。我现在想在每种情况下找到真实分布与估计量定义的分布之间的总变异距离。
由于 30 号和 10 号箱有超过 2 亿种可能的排列方式，我认为使用理论定义不是一个好主意。
包 distrEx has a function "TotalVarDist()", but it can only be used with distributions defined in the distr package, and multinomial is not one of them. There are directions for defining them (see here and here) 但选项是通过显式列出支持来定义离散分布（同样，我认为这不是一个好的选择，因为支持的大小超过 200百万）或使用与创建 distr 包相同的方法从头开始，这超出了我目前的能力。
关于如何使用提到的包或以完全不同的方式执行此操作的任何想法？
我的回答是关于如何使用基数 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 的特定单元格为条件的计数的协方差矩阵时遇到了困难。如果你认为你可以有所作为，我可以分享我的进步。