r - R中连续变量的Tsallis熵

Question

Tsallis entropy对于离散变量定义为:

H[p,q] = 1/(q-1) * (1 - sum(p^q))

连续变量的 Tsallis 熵定义为:

H[p,q] = 1/(q-1) * (1 - int((p(x)^q dx)

哪里 p(x)是数据的概率密度函数， int是不可或缺的。

我正在尝试在 R 中实现 Tsallis 熵。

假设我有以下数据(由 beta 函数生成，但让我们考虑分布未知)

set.seed(567)
mystring <- round(rbeta(500, 2,4), 2)

离散变量的 Tsallis 熵为:

freqs <- table(mystring) / 500
q = 3
H1 <- 1/(q-1) * (1 - sum(freqs^q))
[1] 0.4998426

我现在想计算连续变量的 Tsallis 熵:

PDF <- density(mystring)
library(sfsmisc)
xPDF <- PDF$x
yPDF <- PDF$y
H1 <- 1/(q-1) * (1 - integrate.xy(xPDF, yPDF^q))
[1] -0.6997353

正如我所料，这两个结果是不同的。但为什么如此不同？
我的主要问题是:计算连续变量的 Tsallis 熵的代码是否正确？请记住，我假设分布未知。

Answer 1

所以首先，这是一个统计问题。我鼓励您在 stats.stackexchange.com 上询问它，您可能会得到更好的答案。

话虽如此，您为什么假设值应该相同？您正在从 beta 分布中抽取大小为 n (n = 500) 的随机样本，并尝试通过计算大小为 dx 的 k 个 bin 中的每个 bin 中的观测值的分数(此处为 dx = 0.01 和 k ~ 100)来对其进行离散化。一般来说，每个 bin 中的分数将取决于 k，如

p_i = p_i^o / k

其中 pio 是某些基线 k = ko 的概率向量。换句话说，您拥有的 bin 越多(较小)，每个 bin 的观察次数就越少。您可以通过绘制具有不同 k 的直方图(使用 breaks=k )来看到这一点。

par(mfrow=c(1,3))
hist(mystring,breaks=10,  ylim=c(0,100))
hist(mystring,breaks=50,  ylim=c(0,100))
hist(mystring,breaks=100, ylim=c(0,100))

您的 freqs矢量是 Frequency/500 ，但 k 的效果是一样的。箱数当然等于k，所以

sum( p_i ) = 1

独立于 k。但是在计算 Tsallis 熵时，您不是在对 pi 求和，而是在对 piq 求和(在您的情况下为 q=3)。所以

sum( p_i^q ) ~ sum( [ p_i^o/k ]^q ) ~ (1 / k^q) * sum( [ p_i^o ]^q )

由于您对 k 项求和，当 q = 1 时，结果将不取决于 k，但对于任何其他 q，总和将取决于 k。换句话说，从离散化连续分布计算的 Tsallis 熵将取决于用于离散化的 bin 大小。

为了具体说明这一点，请考虑具有 10 个 bin 的离散化 U[0,1]。这是一个长度为 10 的向量，所有元素 = 0.1。在您的示例中使用 q=3，

k <- 10
p <- rep(1/k,k)
sum(p^q)
# [1] 0.01

现在考虑 100 个 bin 的情况。这里 p 是一个长度为 100 的向量，所有元素都 = 0.01。

k <- 100
p <- rep(1/k,k)
sum(p^q)
# [1] 1e-04

最后考虑连续分布。 U[0,1] = 1 on (0,1), 0 其他地方的 pdf，所以积分是 int(1^3 dx) = 1。

f <- function(x) dunif(x)^q
integrate(f,0,1)$value
# 1

最后，我们可以证明积分您的经验密度函数(基于 dbeta)与直接积分分布函数给出的答案大致相同:

library(sfsmisc)
PDF <- density(mystring)
H2 <- 1/(q-1) * (1 - integrate.xy(PDF$x, PDF$y^q))
H2
# [1] -0.6997353
g <- function(x) dbeta(x,2,4)^q
H3 <- 1/(q-1) * (1 - integrate(g,-Inf,Inf)$value)
H3
# [1] -0.8986014