2018-02-02

エントロピーと幾何平均

$\int p\log{p} dx$ という積分がある。pが確率密度分布であるとき、この積分はKLdの計算の基礎となっているし、 $-\int p\log{p} dx$ はエントロピーである。
$\log{p}$ の期待値とも読めるこの $\int p \log{p} dx$ であるが、離散的な場合の $\sum_{i=1}^n p_i \log{p_i}$ は次のように式変形できる(こちら)
$e^{\sum_{i=1}^n p_i \log{p_i}} = \prod_{i=1}^n p_i^{p_i}$
この右辺が重み付き幾何平均であることが知られており、結局 $\int p \log{p} dx$ は確率変数に関する重み付き幾何平均の対数であることがわかる
以下、重み付き幾何平均について、一般化平均・重み付き一般化平均とともに記載する
一般化平均は $M_q(x)=(\frac{1}{n}\sum_{i=1}^n x_i^q)^{\frac{1}{q}}$ で表される
ただし、 $q=0$ のときには、この式ではだめで、 $M_0(x) = (\prod_{i=1}^n x_i)^{\frac{1}{n}}$ であることが知られている。qについて極限をとった式である
ちなみに、q=0が幾何平均、q=1が算術平均、q=-1が調和平均、q=2がquadratic mean
参考
参考その２
Rで確認を取れば以下のようになる

# generalized mean

my.gen.mean <- function(x,p){
	if(p==0){
		ret <- (prod(x))^(1/length(x))
	}else{
		ret <- ((1/length(x)) * sum(x^p))^(1/p)
	}
	return(ret)
}

ps <- seq(from=-1,to=3,length=100)
ps <- sort(c(ps,0,1))

x <- runif(n)
out1 <- rep(0,length(ps))
for(i in 1:length(out1)){
	out1[i] <- my.gen.mean(x,ps[i])
}

plot(ps,out1)

また、重み付き一般化平均というものもある。値に重みがあり、それは割合(足すと1)であるような場合である
$M_{w,q}(x,w) = (\sum_{i=1}^n w_i x_i^q)^{\frac{1}{q}$
この場合もやはり、q=0のときは極限の式で書いておく必要がある。 $M_{w,q=0}(x,w) = \prod_{i=1}^n x_i^{w_i}$
Rで関数を作れば以下の通り

# weighted generalized mean
w <- runif(n)
w <- w/sum(w)
my.gen.mean2 <- function(x,w,p){
	if(p==0){
		ret <- prod(x^w)
	}else{
		ret <- (sum(w*x^p))^(1/p)
	}
	return(ret)
}

out2 <- rep(0,length(ps))
for(i in 1:length(out2)){
	out2[i] <- my.gen.mean2(x,w,ps[i])
}
plot(ps,out2)

さて。
確率分布を考え、そのばらつきの具合を評価するとき、一般化平均に加えて、次の条件が加わる。 $\sum_{i=1}^n x_i = 1$
また、重み付き一般化平均の場合には、 $x_i=w_i$ である
したがって、確率分布の場合の重み付き一般化平均は
$M_{w,q}(x,w=x)(x) = (\sum_{i=1}^n w_i x_i^q)^{\frac{1}{q}}=(\sum_{i=1}^n x_i^{q+1})^{\frac{1}{q}}$
q=0のときは $M_{w,q=0}(x,w=x)(x) = \prod_{i=1}^n x_i^{w_i}=\prod_{i=1}^n x_i^{x_i}$
重み付き幾何平均の対数が $\sum p \log{p}$ であることを関数で明瞭にすれば

my.plogp <- function(p){
	log(my.gen.mean2(p,p,0))
}
-であり、エントロピーは[tex:H=-\sum p \log(p)]であるので
>|r|
my.entropy <- function(p){
	-my.plogp(p)
}

似たような指標にDiversity Indexというのがある。全体がいくつのカテゴリにどのような割合で分配されているかの指標である
簡単に言うと、重み付き幾何平均の逆数である。したがって $\frac{1}{M_{w,q}(x,w=x)}$ であるのだが
慣習から、 $^uD(x) = \frac{1}{M_{w,u-1}(x,w=x)$ というように、 $u,u-1$ のずれがある
また、Diversity indexのWikiサイトには複数の式が挙げられていて、u,u-1問題もあいまって、わかりにくい。以下のソースでは、Wikiの式表現に沿ってDiversity indexを計算するとともに、それを重み付き幾何平均と対応付けている

n <- 10
p <- runif(n)
#p <- rep(1/n,n)
p <- p/sum(p)
1/prod(p^p)
exp(-sum(p*log(p)))

my.qD <- function(p,q){
	if(q==1){
		ret1 <- 1/prod(p^p)
		ret2 <- exp(-sum(p*log(p)))
	}else{
		ret1 <- 1/((sum(p*p^(q-1)))^(1/(q-1)))
		ret2 <- (sum(p^q))^(1/(1-q))
	}
	
	return(c(ret1,ret2))
}
q <- 1.5
my.qD(p,q)

qs <- seq(from=0,to=2,length=200)
qs <- sort(c(qs,1))
#qs <- qs[-1]
qD.out <- rep(0,length(qs))
for(i in 1:length(qs)){
	qD.out[i] <- my.qD(p,qs[i])[1]
}

plot(qs,qD.out,type="l")

# generalized mean

my.gen.mean <- function(x,p){
	if(p==0){
		ret <- (prod(x))^(1/length(x))
	}else{
		ret <- ((1/length(x)) * sum(x^p))^(1/p)
	}
	return(ret)
}

ps <- seq(from=-1,to=3,length=100)
ps <- sort(c(ps,0,1))

x <- runif(n)
out1 <- rep(0,length(ps))
for(i in 1:length(out1)){
	out1[i] <- my.gen.mean(x,ps[i])
}

plot(ps,out1)


# weighted generalized mean
w <- runif(n)
w <- w/sum(w)
my.gen.mean2 <- function(x,w,p){
	if(p==0){
		ret <- prod(x^w)
	}else{
		ret <- (sum(w*x^p))^(1/p)
	}
	return(ret)
}

out2 <- rep(0,length(ps))
for(i in 1:length(out2)){
	out2[i] <- my.gen.mean2(x,w,ps[i])
}
plot(ps,out2)


my.qD2 <- function(p,q){
	ret <- 1/my.gen.mean2(p,p/sum(p),q-1)
	return(ret)
}

out3 <- rep(0,length(ps))
for(i in 1:length(out3)){
	out3[i] <- my.qD2(x,ps[i])
}
plot(ps,out3)