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

  • こちらの続き
  • 線形回帰と分散分析とで順序関係が等しくなる統計量のことをやっていた
  • 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))))

検定・推定・情報量・エントロピー

  • こちら(ヒストグラムと情報量)の続き
    • これは、ヒストグラムを推定分布とみなして、その分布の情報量のこと、それと分割表の独立を仮定した生起確率・独立性検定統計量のことをメモした記事
  • 最大エントロピー原理(Wiki)の記事にあるように「情報を得る前」の不確かさが「情報を得る」ことによって、「情報を得た後」の不確かさに変化する。不確かさの減分が「情報の量」を表すものとして、「情報量」とする
  • 説明変数と目的変数があって、目的変数がいろいろな値をとるとき、「どんな値をとるか」は「不確か」。そのとき「説明変数」を観測すると、「目的変数」がとる値のいろいろさ加減(「不確かさ」)が変化する。その減分が「説明変数」を観察することの情報量。もし、説明変数を観察すると、目的変数が一意に決まるなら、不確かさは0にまで減じたことになり、この説明変数の観察の情報量は、観察前の目的変数の「不確かさ」そのもの(に相当する情報量)のすべてが、情報として得られたことになる。
  • 2x2分割表の場合
    • 4つの観測度数n_{ij}が得られる。観察に基づく4つの複合カテゴリの生起確率の最尤推定値は観測度数の割合に比例したものとなる(p_{ij}=\frac{n_{ij}}{\sum n_{ij}})。そのような複合4カテゴリの確率分布について-\sum_{i=1}^2 \sum_{j=1}^2 p_{ij}\log{p_{ij}を計算することは、2軸の同時確率分布のエントロピーの計算に相当する
    • 一方、この2x2表について独立性の尤度比検定をすることを考える。対立仮説での対数尤度は\sum_{i=1}^2\sum_{j=1}^2 n_{ij} \log{p_{ij}}である。また帰無仮説での対数尤度は\sum_{i=1}^2\sum_{j=1}^2 n_{ij} \log{\hat{p}_{ij}}であり、尤度比検定では、両者の差の2倍を自由度1のカイ二乗分布に照らして評価するわけであるが、周辺度数が等しい表同士について、統計量の大小のみを、考える場合には、対立仮説のそれは不要になる。それは、n_{ij}の取り方によらず、等しくなるからである。そうすると尤度比検定の統計量の大小順序と4複合カテゴリの推定生起分布のエントロピーの大小順序が等しくなることがわかる

m <- matrix(c(3,4,5,6),2,2)
chisq.test(m,correct=FALSE)

N <- sum(m)
sum(m/N*log(m/N))

n.iter <- 1000
n <- 1000
r <- c(400,600)
c <- c(200,800)
ps <- chisqs <- chisqs2 <- infs <- fps <- x11s <- rep(0,n.iter)
tmp.e <- outer(r,c,"*")/n
tabls <- r2dtable(n.iter,r,c)
for(i in 1:n.iter){
	tmp <- tabls[[i]]
	x11s[i] <- tmp[1,1]
	fout <- fisher.test(tmp)
	fps[i] <- fout$p.value
	chout <- chisq.test(tmp,correct=FALSE)
	ps[i] <- chout$p.value
	chisqs[i] <- chout$statistic
	tmp2 <- tmp[which(tmp!=0)]
	infs[i] <- sum(tmp2/n*log(tmp2/n))
	chisqs2[i] <- 2*n*(infs[i] - sum(tmp2/n*log(tmp.e/n)))
}

plot(chisqs,infs)
plot(fps,infs)

plot(x11s,chisqs)
plot(x11s,infs)

plot(chisqs2,infs)
  • nxm分割表の場合
    • 名義尺度であるなら2x2の場合と同じ説明で、尤度比検定統計量と推定分布のエントロピーの順序が同じであることがわかる
  • 1元配置分散分析の場合
    • 説明変数が2カテゴリ(2群)で目的変数が量的である場合を考える
    • 分散分析は、目的変数の全標本の分散を、群内分散と群間分散に分けて考える
    • 説明変数の情報が得られる前の目的変数のばらつきをエントロピーで表したい。分散の大きさとエントロピーの大きさには結構よい順序関係がある。特に、正規分布の場合には、正規分布標準偏差\sigmaを用いて、そのエントロピー\log{(\sigma\sqrt{2\pi e})}という関係にあるように順序関係が完全に一致している
    • 今、説明変数の情報が得られると、目的変数は、説明変数1に属するか、説明変数2に属するかによって、異なる値の分布を取るだろうと想定されることになる。説明変数のカテゴリ別の分布のエントロピーも観測値値から得られる分散で近似(正規分布近似をすれば、順序は同じ)である。ここで、説明変数のカテゴリ別割合(そのもの、もしくはその推定値)で重みづけをしたエントロピーが説明変数情報を得た後のエントロピーの期待値である。この説明変数によって変化したエントロピーと、説明変数を得る前のエントロピーの差が、説明変数を観察することによって得られた情報量に相当するのだが、全標本分散から説明変数のカテゴリ別の分散をすべて差し引くということは、全標本分散から群内分散を差し引くことであり、その残りの部分は「群間分散」になっている。したがって、説明変数を観察することの情報量は「群間分散」に相当していることになる(ただし、エントロピーでは分布関数の積分になっているのに対して、分散の計算では、総標本数の影響が出ていることは分割表における対数尤度とエントロピーとの関係と同様である

N <- 200
n <- 60
pheno<-c(rep(0,n),rep(1,N-n))
X <- rnorm(N)
n.iter <- 1000
s1 <- s2 <- v1 <- v2 <- rep(0,n.iter)
resid.e <- en1 <- en2 <- rep(0,n.iter)
En <- log(sqrt(var(X))*sqrt(2*pi*exp(1)))
for(i in 1:n.iter){
	tmp.X <- sample(X)
	x <- tmp.X[1:n]
	y <- tmp.X[(n+1):N]
	m.x <- mean(x)
	m.y <- mean(y)
	s1[i] <- cor.test(tmp.X,pheno)[[1]] # 相関
	v1[i] <- var(x)*(n-1)
	v2[i] <- var(y)*(N-n-1)
	s2[i] <- var(X)*(N-1)-(v1[i]+v2[i]) # 群間分散
	en1[i] <- log(sqrt(var(x))*sqrt(2*pi*exp(1)))
	en2[i] <- log(sqrt(var(y))*sqrt(2*pi*exp(1)))
	resid.e[i] <- En-(en1[i]+en2[i])
}
par(mfcol=c(2,2))
plot(s1,s2)
plot(s1^2,s2)
abline(0,1,col=2)

plot(v1,v2)

plot(v1,en1)
par(mfcol=c(1,1))
  • x軸もy軸も量的にするとどうなる?
    • すでに分散分析の統計量と相関係数検定の統計量t = r \sqrt{\frac{df}{1-r^2}}との順序が同じことは上で示した
    • であるから、相関係数検定の統計量(それは相関係数(の2乗)と順序が同じなのだが)についてエントロピーとの関係がわかればよい
    • 相関係数の分子\sum_i (x_i-\bar{x})(y_i-\bar{y})を眺めてみる
    • 分散分析の方でもみたように、各要素と平均との差について、全要素を足し合わせるという作業は、すべての要素の差について足し合わせることと「順序的」に同じなので、\sum_{i,j} (x_i-x_j)(y_i-y_j)と順序的に同じ(あとで計算機的にそうなることを確かめればよい)
    • さらにこの式をじっと見る
    • [tex:\sum_{i,j} *1^2 ]
    • = \sum_{i,j} (x_i-x_j)^2+(y_i-y_j)^2+2(x_i-x_j)(y_i-y_j)が隠れていることがわかる
    • さらに、左辺は[tex:\sum_{i,j}*2^2]であって、これは、2点がx+y=Cという直線のどこに乗るかで、その(C_i-C_j)^2という値のことになる
    • また、右辺は、xに関する分散のようなものとyに関する分散のようなものの和と、今、着目している\sum_{i,j} (x_i-x_j)(y_i-y_j)になっている。ここでx,yのそれぞれの分散のようなものは、x,yについて「すでに与えられている」状況では変化がないから、\sum_{i,j} (x_i-x_j)(y_i-y_j)(C_i-C_j)^2との順序は同じになることがわかる
    • これは結局、2点の距離の2乗について、相関を考えている方向(ベクトル(1,1)の方向成分と、それと直交する成分とに分けることで、相関を考えている成分の多寡を評価していることを意味している
    • 2次元の場合は、負の相関(ベクトル(1,-1))がもう一つの方向としてあり(ただし、これは計算しなくても正の相関の裏返しとして算出される)し、次元が上がったら(1,1,1),(1,1,-1),(1,-1,-1),...のように複数の考慮するべき方向が出てくることも意味している
    • ここでは、平均・分散につながる2次の計算しかしていないので、2次以外の情報はゼロ、従って、2次元正規分布を想定していることになる。その2次元正規分布の分散は直交する2軸の分散の和となるが、そのような2軸は楕円の長径と短径の方向であるときにそのようになる。与えられたx軸、y軸の分散の和が観察している楕円型の分散に等しくならないときには、その差の分だけ、「情報が減じて」いるその量の多寡を使って検定していることになる(らしい)
    • \sum_{i,j} (x_i-x_j)(y_i-y_j):この式が群間分散と同様のものであるのはなぜかというと、説明変数Xについてもし同じカテゴリに属していればx_i=x_jであるから、(x_i-x_j)(y_i-y_j)が0になって「加えない」のに対し、x_i \ne x_jである場合には、何かしらx_i-x_jという値の重みでy_i-y_jが足されることになる。ただし、説明変数Xのすべてのカテゴリは相互に対称的なのでx_i-x_ji,jによらず(順番にもよらず)同じ値を持つものと仮定するから(それに1というありきたりな値を仮定するのと変わらない)、そういう値になっている。y_i-y_jが「方向付き距離」だから、これはz検定対象(正規分布で考える)…「距離の2乗」にすればカイ二乗分布で考える…と言ったところだろうか
    • コメントも書き込まないけれど、これを確かめるRのソース
sum.sq <- function(x){
	sum((outer(x,x,"-"))^2)/length(x)
}

sum.sq.d <- function(x){
	sum(dist(x,method="manhattan")^2)
}

tmp.sum <- 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(x[i,]-x[j,]))^2
		}
	}
	sum(m)
}
tmp.sum.2 <- function(x){
	if(!is.matrix(x)){
		x <- matrix(x,ncol=1)
	}
	a.x <- apply(x,1,sum)
	sum((outer(a.x,a.x,"-"))^2)
}

N <- 200
n <- 60
X <- rnorm(N)

pheno<-c(rep(0,n),rep(0.00001,N-n))
pheno <- runif(N)
n.iter <- 100
s1 <- s2 <- s3 <- s4 <- s5 <- v1 <- v2 <- rep(0,n.iter)
resid.e <- en1 <- en2 <- rep(0,n.iter)
En <- log(sqrt(var(X))*sqrt(2*pi*exp(1)))
for(i in 1:n.iter){
	tmp.X <- sample(X)
	x <- tmp.X[1:n]
	y <- tmp.X[(n+1):N]
	m.x <- mean(x)
	m.y <- mean(y)
	s1[i] <- cor.test(tmp.X,pheno)[[1]] # 相関
	s3[i] <- cor.test(tmp.X,pheno)[[4]] # 相関
	#s4[i] <- sum.sq.d(cbind(tmp.X,pheno))-(sum.sq.d(tmp.X)+sum.sq.d(pheno))
	s4[i] <- tmp.sum.2(cbind(tmp.X,pheno))-(tmp.sum.2(tmp.X)+tmp.sum.2(pheno))
	s5[i] <- tmp.sum.2(cbind(tmp.X,pheno))
	#v1[i] <- var(x)*(n-1)
	#v2[i] <- var(y)*(N-n-1)
	#s2[i] <- var(X)*(N-1)-(v1[i]+v2[i]) # 群間分散
	v1[i] <- sum.sq(x)
	v2[i] <- sum.sq(y)
	s2[i] <- sum.sq(X)-(v1[i]+v2[i]) # 群間分散
	#s2[i] <- (v1[i]+v2[i]) # 群間分散
	en1[i] <- log(sqrt(var(x))*sqrt(2*pi*exp(1)))
	en2[i] <- log(sqrt(var(y))*sqrt(2*pi*exp(1)))
	resid.e[i] <- En-(en1[i]+en2[i])
}
par(mfcol=c(2,2))
plot(s1,s2)
plot(s1^2,s2)
plot(s3,s2)
plot(s1,s4)
#abline(0,1,col=2)

#plot(v1,v2)

#plot(v1,en1)
par(mfcol=c(1,1))
  • 3群以上のone-way ANOVA
    • 群を独立な説明変数として、(群数-1)次元平面上に目的変数の雲があって、そkに「平面」を線形回帰して検定に持ち込む
    • Rではlm()関数で線形回帰をして、その推定係数を用いて、anova()関数で検定統計量(F統計量)に照らして検定する
    • 目的変数の分散と説明ファクタごとの分散(群内分散)とに分解することを考えれば、線形回帰経由の検定統計量と群内分散の和とが等しい順序になる
    • Rでやろう

sum.sq <- function(x){
	sum((outer(x,x,"-"))^2)/length(x)
}


ns <- sample(10:50,5)
N <- sum(ns)
X <- rnorm(N)
Pheno <- matrix(0,N,length(ns))
Pheno.factor <- c()
cnt <- 1
for(i in 1:length(ns)){
	Pheno[cnt:(cnt+ns[i]-1),i] <- 1
	cnt <- cnt+ns[i]
	tmp <- paste("",i)
	Pheno.factor <- c(Pheno.factor,rep(tmp,ns[i]))
}

n.iter <- 100

a <- b <- rep(0,n.iter)
for(i in 1:n.iter){
	tmp.X <- sample(X)
	tmp.a <- lm(tmp.X ~ Pheno.factor)
	tmp.a2 <- anova(tmp.a)
	a[i] <- tmp.a2[[4]][1]
	S.all <- sum.sq(tmp.X)
	S.intra <- rep(0,length(ns))
	cnt <- 1
	for(j in 1:length(ns)){
		S.intra[j] <- sum.sq(tmp.X[cnt:(cnt+ns[j]-1)])
		cnt <- cnt+ns[j]
	}
	b[i] <- sum(S.intra)
}

plot(a,b)

*1:x_i-x_j)+(y_i-y_j

*2:x_i+y_i) - (x_j+y_j