2013-06-19

エントロピーで関連

昨日のさらに続き
２軸でいわゆる相関係数を取る場合と、エントロピーとを見る場合とで、どんな違いが出るかを、

にあるような例でやってみよう

この絵はWikipediaのCorrelationの例の図
- 上段は傾き45度でばらつきを変えている
- 中段はばらつきは同じで傾きを変えている
- 下段はパターンはあるものの、いわゆる相関がとれるとは限らない・とれないもの
エントロピーベースで評価すると、ばらつきの強弱を反映する、傾きはまったく関係ない、視覚パターンもそれなりに評価する

library(FNN)
Entropy.perm <- function(X,n.perm=100,k=1){
	Ent.ori <- entropy(X,k=k)[k]
	Ent.perm <- rep(0,n.perm)
	tmp <- cor.test(X[,1],X[,2])
	Cor.ori <- tmp[[1]]
	#cor.p.ori <- tmp[[3]]
	Cor.perm <- Ent.perm
	for(i in 1:n.perm){
		tmp.X <- X
		tmp.X[,1] <- sample(tmp.X[,1])
		Ent.perm[i] <- entropy(tmp.X,k=k)[k]
		tmp <- cor.test(tmp.X[,1],tmp.X[,2])
		Cor.perm[i] <- tmp[[1]]
	}
	p <- length(which(Ent.perm < Ent.ori))/n.perm
	cor.p <- length(which(Cor.perm < Cor.ori))/n.perm
	
	return(list(Ent.ori=Ent.ori,Ent.perm=Ent.perm,p.value=p,cor.p=cor.p,Cor.perm=Cor.perm,Cor.ori=Cor.ori))
}
n <- 100

X <- matrix(rnorm(n*2),ncol=2)

ys <- seq(from=0.3,to=0.8,length=20)
theta <- pi/4
#par(mfcol=c(4,length(ys)/4))
Rot <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)
cors <- KLs <- Ents <- Cors <- rep(0,length(ys))
#KLs.list <- list()
Ent.list <- list()
Ent.p <- cor.perm.p <- cor.p <- Ents
sorted.X <- cbind(sort(X[,1]),sort(X[,2]))
for(i in 1:length(ys)){
	tmp.X <- X
	tmp.X[,2] <- tmp.X[,2]*ys[i]
	tmp.X <- tmp.X %*% t(Rot)
	#plot(tmp.X)
	cor.out <- cor.test(tmp.X[,1],tmp.X[,2])
	cors[i] <- cor.out[[1]]
	cor.p[i] <- cor.out[[3]]
	#KLs[i] <- KL.dist(tmp.X[,1],tmp.X[,2],k=1)[1]
	#KLs.list[[i]] <- KL.dist.perm(sorted.X,tmp.X,n.perm=100)
	#KLs[i] <- KLs.list[[i]]$p
	Ent.list[[i]] <- Entropy.perm(tmp.X,k=1)
	Ents[i] <- Ent.list[[i]]$Ent.ori
	Cors[i] <- Ent.list[[i]]$Cor.ori
	Ent.p[i] <- Ent.list[[i]]$p.value
	cor.perm.p[i] <- Ent.list[[i]]$cor.p
	
}
#par(mfcol=c(1,1))
dev.new()
par(mfcol=c(2,2))
plot(cors,Ents)
plot(cor.p,Ent.p)
plot(cor.p,cor.perm.p)
plot(cor.perm.p,Ent.p)
par(mfcol=c(1,1))



y <- 0.2
thetas <- seq(from=-pi/4,to=pi/4,length=20)
#par(mfcol=c(4,length(ys)/4))
Rot <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)
cors <- KLs <- Ents <- Cors <- rep(0,length(thetas))
#KLs.list <- list()
Ent.list <- list()
Ent.p <- cor.perm.p <- cor.p <- Ents
sorted.X <- cbind(sort(X[,1]),sort(X[,2]))
for(i in 1:length(thetas)){
	theta <- thetas[i]
	tmp.X <- X
	tmp.X[,2] <- tmp.X[,2]*y
	Rot <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)

	tmp.X <- tmp.X %*% t(Rot)
	#plot(tmp.X)
	cor.out <- cor.test(tmp.X[,1],tmp.X[,2])
	cors[i] <- cor.out[[1]]
	cor.p[i] <- cor.out[[3]]
	#KLs[i] <- KL.dist(tmp.X[,1],tmp.X[,2],k=1)[1]
	#KLs.list[[i]] <- KL.dist.perm(sorted.X,tmp.X,n.perm=100)
	#KLs[i] <- KLs.list[[i]]$p
	Ent.list[[i]] <- Entropy.perm(tmp.X,k=1)
	Ents[i] <- Ent.list[[i]]$Ent.ori
	Cors[i] <- Ent.list[[i]]$Cor.ori
	Ent.p[i] <- Ent.list[[i]]$p.value
	cor.perm.p[i] <- Ent.list[[i]]$cor.p
	
}
#par(mfcol=c(1,1))
#dev.new()
#par(mfcol=c(2,2))
plot(cors,Ents)
#plot(cor.p,Ent.p)
#plot(cor.p,cor.perm.p)
#plot(cor.perm.p,Ent.p)
par(mfcol=c(1,1))

n <- 1000
ss <- seq(from=0,to=2,length=100)
ss <- ss[-1]
x <- rnorm(n)
ent.ps <- cor.ps <- rep(0,length(ss))
for(i in 1:length(ss)){
	y <- cos(x*pi/2) + rnorm(length(x),0,ss[i])
	plot(x,y)
	tmp <- Entropy.perm(cbind(x,y),n.perm=100)
	ent.ps[i] <- tmp$p.value
	cor.ps[i] <- tmp$cor.p
}
matplot(cbind(ent.ps,cor.ps),type="l")

n <- 100
t <- runif(n)*2*pi
x <- cos(t)
y <- sin(t)
ss <- seq(from=0,to=0.2,length=100)
ss <- ss[-1]
ent.ps <- cor.ps <- rep(0,length(ss))
for(i in 1:length(ss)){
	x <- x + rnorm(length(x),0,ss[i])
	y <- y + rnorm(length(x),0,ss[i])
	plot(x,y)
	tmp <- Entropy.perm(cbind(x,y),n.perm=100)
	ent.ps[i] <- tmp$p.value
	cor.ps[i] <- tmp$cor.p
}
matplot(cbind(ent.ps,cor.ps),type="l")

n <- 1000
x <- runif(n)
y <- runif(n)
X <- cbind(x,y)
thetas <-seq(from=0,to=pi/4,length=10)
ent.ps <- cor.ps <- rep(0,length(thetas))
for(i in 1:length(thetas)){
	Rot <- matrix(c(cos(thetas[i]),sin(thetas[i]),-sin(thetas[i]),cos(thetas[i])),2,2)
	tmp.X <- X %*% Rot
	plot(tmp.X)
	tmp <- Entropy.perm(tmp.X,n.perm=100)
	ent.ps[i] <- tmp$p.value
	cor.ps[i] <- tmp$cor.p
}
matplot(cbind(ent.ps,cor.ps),type="l")