2013-06-15

検定・推定・情報量・エントロピー(続き)

こちらの続き
線形回帰と分散分析とで順序関係が等しくなる統計量のことをやっていた
１元配置分散分析で２群のときは、線形回帰そのものだった
１元配置分散分析で３以上の群のときは群内分散・群間分散の割り振り具合と群の数(もしくはそれより１つ少ない群の数)を説明変数とした線形回帰のあてはまりで考えればよいのだった
その場合に、いわゆる普通の相関係数のようなものがあるだろうか
- ２群の場合には、説明変数が１次元空間の２方向
- n群の場合には説明変数がn-1次元空間のn方向であると考えると、「相関係数」はベクトル様になる
- ２群の場合の相関係数が「スカラー」なのは、１次元空間の２方向で表されたベクトルはそのどちらか片方の大きさを問題にすればよいから
- n群の場合の相関係数ベクトルについては、自由度をn-1で考えるならば、この相関係数ベクトルの長さの２乗の和がn群の自由度n-1での相関係数として使えそうだ
この量と「情報量」の関係がわからないとこの記事としては意味がないのだが…
- 説明変数と目的変数が作る分布のばらつきの大きさと説明変数のそれと目的変数のそれとの和との差が情報量
- 平均と分散しか気にしないとき〜正規分布を仮定しているとき〜結果として正規分布を仮定しているわけでそれしか知らないからそうやっているというナイーブなものかもしれない〜ときには、ペアワイズで目的変数の差に「応じて(向きも含めて)」説明変数の差(向きも含めて)を重みづけることで、方向ありのばらつきベクトルが標本ペアの数だけできる。その総和は、説明変数と目的変数の関係の強さと関係した値になる
- どの向きに両説明変数が関係しているかに興味がないなら、その総和は向きを無視して大きさだけで足し合わせるのがよかろう、と、そういう説明になるのだろう
そうすると、目的変数が「ベクトルであるとき」には「ベクトルの違い」を「説明変数の違い」に応じて足し合わせればよく、自由度がフルなら、ノルムの足し合わせになるか…
また、説明変数の違いを線形に考えるか、そうしないかによって、さらに一般化できるだろう

sum.sq <- function(x){
	sum((outer(x,x,"-"))^2)/length(x)
}
tmp.sum.3 <- function(x){
	if(!is.matrix(x)){
		x <- matrix(x,ncol=1)
	}
	m <- matrix(0,length(x[,1]),length(x[,1]))
	for(i in 1:(length(x[,1])-1)){
		for(j in 2:length(x[,1])){
			m[i,j] <- m[j,i] <- (sum(abs(x[i,]-x[j,])))^2
		}
	}
	sum(m)
}
CategoryVector<-function (nc = 3){
	df <- nc - 1
	d <- df + 1
	diagval <- 1:d
	diagval <- sqrt((df + 1)/df) * sqrt((df - diagval + 1)/(df - diagval + 2))
	others <- -diagval/(df - (0:(d - 1)))
	m <- matrix(rep(others, df + 1), nrow = df + 1, byrow = TRUE)
	diag(m) <- diagval
	m[upper.tri(m)] <- 0
	as.matrix(m[, 1:df])
}

make.simplex <- function(k){
	cv <- CategoryVector(k)
	rbind(t(cv*sqrt(1-1/k)),rep(1/sqrt(k),k))
}

make.simplex.multi <- function(r){
	X <- make.simplex(r[1])
	k <- length(r)
	if(k > 1){
		for(i in 2:k){
			X <- X %x% make.simplex(r[i])
		}
	}
	X
}
sum.sq <- function(x){
	sum((outer(x,x,"-"))^2)/length(x)
}

library(MCMCpack)
N <- 100
d <- 3
cv <- CategoryVector(d)
preX.cat <- matrix(0,N,d)
preX.cat[,1] <- 1
preX.cat <- t(apply(preX.cat,1,sample))
X.cat <- preX.cat%*%cv
#preX <- rdirichlet(N,rep(0.01,d))
#X <- preX%*%cv
preX <- preX.cat + rnorm(length(preX.cat),0,0.01)
X <- preX%*%cv
plot(X.cat)
plot(X)

Y <- rnorm(N)

n.iter <- 100
s1 <- s2 <- s3 <- s4 <- s5 <- rep(0,n.iter)
v1 <- v3 <- matrix(0,n.iter,d)
v2 <- v4 <- matrix(0,n.iter,d-1)
v1.2 <- v3.2 <- v1
v2.2 <- v4.2 <- v2
dst1 <- array(0,c(N,N,d))
dst2 <- array(0,c(N,N,d-1))
dst3 <- array(0,c(N,N,d))
dst4 <- array(0,c(N,N,d-1))

for(i in 1:N){
	for(j in 1:N){
		dst1[i,j,] <- preX.cat[i,]-preX.cat[j,]
		dst2[i,j,] <- X.cat[i,]-X.cat[j,]
		dst3[i,j,] <- preX[i,]-preX[j,]
		dst4[i,j,] <- X[i,]-X[j,]
	}
}

for(i in 1:n.iter){
	shY<-sample(Y)
	tmp1 <- rep(0,d)
	for(j in 1:d){
		ss <- which(preX.cat[,j]==1)
		tmp1[j] <- sum.sq(shY[ss])
	}
	s1[i] <- sum(tmp1)
	tmp.a <- lm(shY ~ preX.cat)
	tmp.a2 <- anova(tmp.a)
	s2[i] <- tmp.a2[[4]][1]
	tmp.b <- lm(shY ~ X.cat)
	tmp.b2 <- anova(tmp.b)
	s3[i] <- tmp.b2[[4]][1]
	#dstY <- as.matrix(dist(shY))
	dstY <- (outer(shY,shY,"-"))
	
	for(j in 1:length(v1[i,])){
		v1[i,j] <- sum(dst1[,,j]*dstY)
		v3[i,j] <- sum(dst3[,,j]*dstY)
		if(j<length(v1[i,])){
			v2[i,j] <- sum(dst2[,,j]*dstY)
			v4[i,j] <- sum(dst4[,,j]*dstY)

		}
	}
	#s4[i] <- sum((dst1*dstY))
}
plot(as.data.frame(cbind(s1,s2,s3,apply(v1^2,1,sum),apply(v2^2,1,sum),apply(v3^2,1,sum),apply(v4^2,1,sum))))