- こちらの続き
- 分散がエントロピーと関係があることや、群内分散と群間分散とへの分解が分散分析であること、分散をエントロピーの代わりとすることは正規分布を仮定していることなどを書いてきた
- 正規分布とは違う分布として、多峰性の分布を取り上げる。多峰性分布では分散がエントロピーの代用にならないから
- こんな分布を考える。0-1の範囲の正弦関数を確率密度分布とする。２つの分布は位相がずれている
t <- seq(from=0,to=1,length=100)
k <- 4
theta <- pi/2
Xa <- (sin(t*k*pi)+1)/2
Xb <- (sin(t*k*pi+theta)+1)/2
matplot(cbind(Xa,Xb),type="l")
- こんな確率密度分布からの標本について
- ２群の分散分析
- Nearest-Neighbor KL-distance estimation
- の２つの統計量でパーミュテーショナルなp値を出してみる
- 正弦関数の周期を変数として変える。分散分析は(少なくとも片方の)分布の峰が１つなら違いをp値に反映できるが、峰の数が複数になるとすぐにだめになるのに対して、KL-distの方は結構(標本数が分布をそれなりに体現する程度に多いならば)検出できるようである
- 位相のずれの大きさについて検討してもよいだろうし、標本数とp値の関係について検討してもよいだろう
- KL-distanceの推定には最近接点を用いる方法と近接する複数の点を用いる方法があるが、以下では最近接点を用いている(近接するいくつかの点にする場合には、それをいくつにするかを定める必要も出る)
- k-Nearest Neighbor法の基礎には,k-Nearest Neighbor密度推定がある→こちら
- 横軸にsin(t*k*pi)のkを、縦軸にp値の平均を取ると：
library(FNN)
ks <- c(0,2^(seq(from=-4,to=4,length=4)))
ks <-seq(from=0,to=4,length=10)
ks <- 2^(1:10)
thetas <- seq(from=0,to=1,length=9)*pi/2
thetas <- pi/2
k.thetas <- expand.grid(ks,thetas)
na <- 100
nb <- 100
n.rep <- 100
kl.k <- 1
S.ps <- K.ps <- matrix(0,n.rep,length(ks))
for(rr in 1:n.rep){
S.ori <- rep(0,length(k.thetas[,1]))
K.ori <- matrix(0,length(S.ori),10)
S.p <- S.ori
K.p <- S.p
for(ii in 1:length(k.thetas[,1])){
k <- k.thetas[ii,1]
theta <- k.thetas[ii,2]
xa <- xb <- c()
for(i in 1:na){
loop <- TRUE
while(loop){
r <- runif(2)
if(r[1] < (sin(r[2]*k*pi)+1)/2){
xa <- c(xa,r[2])
loop <- FALSE
}
}
}
for(i in 1:nb){
loop <- TRUE
while(loop){
r <- runif(2)
if(r[1] < (sin(r[2]*k*pi+theta)+1)/2){
xb <- c(xb,r[2])
loop <- FALSE
}
}
}
S.ori[ii] <- sum.sq(xa)+sum.sq(xb)
K.ori[ii,] <- KL.dist(xa,xb,k=kl.k)
n.iter <- 100
S.perm <- rep(0,n.iter)
K.perm <- matrix(0,n.iter,length(K))
for(i in 1:n.iter){
sh <- sample(c(xa,xb))
tmp.xa <- sh[1:na]
tmp.xb <- sh[(na+1):(na+nb)]
S.perm[i] <- sum.sq(tmp.xa)+sum.sq(tmp.xb)
K.perm[i,] <- KL.dist(tmp.xa,tmp.xb,k=kl.k)
}
S.p[ii] <- length(which(S.ori[ii] > S.perm))/n.iter
K.p[ii] <- length(which(K.ori[ii,1] < K.perm[,1]))/n.iter
}
matplot(cbind(S.p,K.p),type="l")
S.ps[rr,] <- S.p
K.ps[rr,] <- K.p
}
S.mean <- apply(S.ps,2,mean)
K.mean <- apply(K.ps,2,mean)
matplot(cbind(S.mean,K.mean),type="l")
boxplot(S.ps)
boxplot(K.ps)
