R中连续变量的Tsallis熵

Question

离散变量的Tsallis熵定义如下：

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熵代码是否正确？请记住，我假设分布是未知的。

- Tommaso

Answer 1

首先，这是一个统计学问题。我建议你在stats.stackexchange.com上提问，那里你可能会得到更好的答案。

话虽如此，为什么你假设这些值应该相同呢？你从Beta分布中随机抽取样本大小n（n=500），并尝试通过计算每个k个大小为dx（这里dx=0.01，k〜100）的存储桶中观测值的比例来离散化它。一般来说，每个桶中的比例将取决于k，如下所示：

p_i = p_i^o / k

其中，p_i^o是某个基线k=k_o的概率向量。换句话说，你拥有的（更小的）桶越多，每个桶中的观测次数就越少。使用breaks=k绘制具有不同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熵时，你不是求和p_i，而是求和p_i^q（在你的例子中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熵将取决于用于离散化的箱体大小。

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

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

现在考虑使用100个箱子的情况。这里p是一个长度为100的向量，所有元素都等于0.01。

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

最后考虑连续分布。U[0,1]的概率密度函数在(0,1)为1，在其他地方为0，所以积分结果为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

- jlhoward