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

  • こちらの続き
  • 線形回帰と分散分析とで順序関係が等しくなる統計量のことをやっていた
  • 1元配置分散分析で2群のときは、線形回帰そのものだった
  • 1元配置分散分析で3以上の群のときは群内分散・群間分散の割り振り具合と群の数(もしくはそれより1つ少ない群の数)を説明変数とした線形回帰のあてはまりで考えればよいのだった
  • その場合に、いわゆる普通の相関係数のようなものがあるだろうか
    • 2群の場合には、説明変数が1次元空間の2方向
    • n群の場合には説明変数がn-1次元空間のn方向であると考えると、「相関係数」はベクトル様になる
    • 2群の場合の相関係数が「スカラー」なのは、1次元空間の2方向で表されたベクトルはそのどちらか片方の大きさを問題にすればよいから
    • n群の場合の相関係数ベクトルについては、自由度をn-1で考えるならば、この相関係数ベクトルの長さの2乗の和が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))))