検定・推定・情報量・エントロピー - ryamadaの遺伝学・遺伝統計学メモ

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

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)

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

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))

ｘ軸もｙ軸も量的にするとどうなる？
- すでに分散分析の統計量と相関係数検定の統計量 $t = r \sqrt{\frac{df}{1-r^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]であって、これは、２点が $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$ との順序は同じになることがわかる
- これは結局、２点の距離の２乗について、相関を考えている方向(ベクトル $(1,1)$ の方向成分と、それと直交する成分とに分けることで、相関を考えている成分の多寡を評価していることを意味している
- ２次元の場合は、負の相関(ベクトル $(1,-1)$ )がもう一つの方向としてあり(ただし、これは計算しなくても正の相関の裏返しとして算出される)し、次元が上がったら $(1,1,1),(1,1,-1),(1,-1,-1),...$ のように複数の考慮するべき方向が出てくることも意味している
- ここでは、平均・分散につながる２次の計算しかしていないので、２次以外の情報はゼロ、従って、２次元正規分布を想定していることになる。その２次元正規分布の分散は直交する２軸の分散の和となるが、そのような２軸は楕円の長径と短径の方向であるときにそのようになる。与えられたｘ軸、ｙ軸の分散の和が観察している楕円型の分散に等しくならないときには、その差の分だけ、「情報が減じて」いるその量の多寡を使って検定していることになる(らしい)
- $\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_j$ は $i,j$ によらず(順番にもよらず)同じ値を持つものと仮定するから(それに１というありきたりな値を仮定するのと変わらない)、そういう値になっている。 $y_i-y_j$ が「方向付き距離」だから、これはz検定対象(正規分布で考える)…「距離の２乗」にすればカイ二乗分布で考える…と言ったところだろうか
- コメントも書き込まないけれど、これを確かめる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))

３群以上のone-way ANOVA
- 群を独立な説明変数として、(群数-１)次元平面上に目的変数の雲があって、そｋに「平面」を線形回帰して検定に持ち込む
- Rではlm()関数で線形回帰をして、その推定係数を用いて、anova()関数で検定統計量(Ｆ統計量)に照らして検定する
- 目的変数の分散と説明ファクタごとの分散(群内分散)とに分解することを考えれば、線形回帰経由の検定統計量と群内分散の和とが等しい順序になる
- 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