2013-09-07

分割表として一般化

決断診断ベイジアンネットワークベータ分布ロジスティック回帰射影多様体微分幾何カルマンフィルタ

昨日の話題の一般化
昨日は2x2表からベータ分布を使って、事後生起確率分布を算出する話と、それを数値計算的に出す話を書いた
数値計算にするにしても、次元が膨大になると大変そう…と思ったり、要因の数が大きくなりすぎると、overfittingがどういう形で効いてくるだろうか、とか、気になることが次から次にでてくるし
また、連続説明変数からのロジスティック回帰への転用とか
従属変数が連続のときのノンパラ・パラのこととか
も気になるわけだけれど、いったん、そのことを忘れて、ベータ分布の活用という観点だけをとにかく、一般化して書くとどうなるか、ということにしぼってみよう

my.marginal.apply <- function(A,v,fun="sum"){
	if(length(v)==0){
		return(fun(A))
	}else{
		return(apply(A,v,fun))
	}
}

my.multi.cat <- function(D,P){
	n <- length(D[1,])
	tmp <- apply(t(D)*2^(0:(n-1)),2,sum)+1
	tmp.1 <- tmp[which(P==0)]
	tmp.2 <- tmp[which(P==1)]
	tab.1 <- tabulate(tmp.1,2^n)
	tab.2 <- tabulate(tmp.2,2^n)
	return(list(A1=tab.1,A2=tab.2))
}
my.p.seq <- function(p.len){
	tmp <- seq(from=0,to=1,length=p.len)
	tmp <- tmp[-1]
	tmp <- tmp[-length(tmp)]
	return(tmp)
}
my.d.beta <- function(D,P,phenotype,p.len=100,pre.p=rep(1/(2^length(D[1,])),2^length(D[1,])),model="all"){
	# 説明変数の数
	n <- length(D[1,])
	# 2^nのセルのパターン
	#X <- expand.grid(rep(list(0:1),n))
	
	if(model=="full"){
		X <- rbind(rep(0,n),rep(1,n))
	}else if(model=="all"){
		X <- expand.grid(rep(list(0:1),n))
	}else if(model=="alternative"){
		X <- matrix(rep(1,n),nrow=1)
	}
	
	As <- my.multi.cat(D,P)
	A1 <- array(As[[1]],rep(2,n))
	A2 <- array(As[[2]],rep(2,n))
	# 確率・尤度を計算するための「成功確率の値」
	p <- my.p.seq(p.len)

	# 説明変数の取り方パターン(2^n通りある)ごとに、その「とり方パターンが真である」という仮説をどれくらいしんじたらよいか、その仮説の場合に特定のphenotypeでの成否確率分布を値ベクトルpについて計算する
	# 説明変数の採用パターンごとの尤度の変化の比(すべての説明変数で2^n通りに分けた場合を基準とする)
	K <- rep(0,length(X[,1]))
	prs.all <- list()
	for(i in 1:length(phenotype[,1])){
		prs.all[[i]] <- matrix(0,length(X[,1]),length(p))
	}
	# 説明変数の採用パターンごとの、特定phenotypeでの成否確率分布
	#prs <- matrix(0,length(K),length(p))
	for(i in 1:length(X[,1])){
		# この仮説で問題にする説明変数を取り出す
		s <- which(X[i,]==1)
		# apply()関数で部分分割表化する
		tmp1 <- my.marginal.apply(A1,s,sum)
		tmp2 <- my.marginal.apply(A2,s,sum)
		# 成否ごとのカウントとその和とが必要
		marg <- tmp1+tmp2
		# 限定した仮説に関して、フルモデルのもとでの尤度にもとづく「比」を計算する
		K[i] <- sum(lgamma(marg+2)-lgamma(tmp1+1)-lgamma(tmp2+1))
		# 仮説においては、特定のセルにおける成否カウントとのその和に基づくベータ分布が必要になる
		# apply()関数で作った部分分割表のどのセルがその特定のセルかを与える
		locs <- rep(1,length(prs.all))
		for(j in 1:length(locs)){
			if(length(s)>0){
				locs[j] <- sum(phenotype[j,s]*2^(0:(length(s)-1)))+1
			}
			prs.all[[j]][i,] <- dbeta(p,tmp1[locs[j]]+1,tmp2[locs[j]]+1)
		}
		#loc <- 1
		#if(length(s)>0){
		#	loc <- sum(phenotype[s]*2^(0:(length(s)-1)))+1
		#}
		#prs[i,] <- dbeta(p,tmp1[loc]+1,tmp2[loc]+1)
	}
	# フルモデルを基準に比を取る
	# 対数で計算していたから、加減計算して、指数関数にかける
	K.2 <- -K+K[length(K)]
	# 仮説の事前確率はすべて同一だったとする
	pre.p <- rep(1/length(X[,1]),length(X[,1]))
	# これが仮説別の事後確率
	post.p <- pre.p * exp(K.2)
	post.p <- post.p/sum(post.p)
	# 仮説別の事後確率で重みづけをして、統合的成功事後確率分布を求める
	PR.all <- matrix(0,length(prs.all),length(p))
	for(j in 1:length(phenotype[,1])){
		PR <- rep(0,length(p))
		for(i in 1:length(pre.p)){
			PR <- PR + post.p[i]*prs.all[[j]][i,]
		}
		PR.all[j,] <- PR
	}

	return(list(As=As,phenotype=phenotype,p=p,pre.p=pre.p,post.p=post.p,prs=prs.all,LL=K.2,est.pr=PR.all))
}
my.d.beta.models <- function(D,P,phenotype,p.len=100,pre.p=rep(1/(2^length(D[1,])),2^length(D[1,])),models=c("all","full","alt")){
	ret <- list()
	for(i in 1:length(models)){
		ret[[i]] <- my.d.beta(D,P,phenotype,p.len=p.len,pre.p=pre.p,model=models[i])
	}
	return(list(models=models,result = ret))
}
n<-4




n.sample <- 1000
D <- matrix(sample(0:1,n*n.sample,replace=TRUE),ncol=n)
P <- sample(0:1,n.sample,replace=TRUE)
phenotype <- expand.grid(rep(list(0:1),n))
p.len <- 100
model="alternative"

tmp.out <- my.d.beta(D,P,phenotype,model=model)

matplot(tmp.out$p,t(tmp.out$est.pr),type="l")
for(i in 1:length(phenotype[,1])){
	par(mfcol=c(1,2))
	matplot(tmp.out$p,cbind(tmp.out$est.pr[i,],t(tmp.out$prs[[i]])),type="l",ylim=c(0,max(tmp.out$prs[[i]])),lwd=c(5,rep(1,length(tmp.out$prs))))
	par(ask=FALSE)
	plot(tmp.out$p,tmp.out$est.pr[i,],type="l",ylim=c(0,max(tmp.out$prs[[i]])))
	
	par(mfcol=c(1,1))
	par(ask=TRUE)
}

library(MCMCpack)
# 説明変数の数
n <- 5
# 2^nのセルのパターン
X <- expand.grid(rep(list(0:1),n))
X
# 2^nパターンごとに、成否の２観察をしたとする
# ここは適当
y <- round(rdirichlet(1,rep(0.3,length(X[,1])*2))*500)
y
#y<- rep(sample(0:2,length(X[,1]),replace=TRUE),2)
y <- y + sample(0:3,length(y),replace=TRUE)
y
# ２列(成否)の行列にする
# 行は説明変数のパターンに対応
Y <- matrix(y,ncol=2)
# 確率・尤度を計算するための「成功確率の値」
p <- seq(from=0,to=1,length=100)
# 確率として0,1を外す
p <- p[-length(p)]
p <- p[-1]
# 長さ2^nのベクトルをn次元のアレイに納め
# apply()関数を使って、説明変数の採用の仕方に関する2^n通りのすべてに関して周辺度数化(問題にする説明変数と問題にしない説明変数に分けて、問題にする説明変数だけに関する部分分割表にする)するための関数
# apply()関数をarrayにかけるとき、『全要素の和』にもたいおうするようにちょっと細工した関数
my.marginal.apply <- function(A,v,fun="sum"){
	if(length(v)==0){
		return(fun(A))
	}else{
		return(apply(A,v,fun))
	}
}
# 成否確率の推定をしたい特定の説明変数パターンを一つ決める
phenotype <- sample(0:1,n,replace=TRUE)

# 説明変数パターン別のカウントデータであるYを成否ごとに異なるn次元アレイに納める
A1 <- array(Y[,1],rep(2,n))
A2 <- array(Y[,2],rep(2,n))

A1
A2

# 説明変数の取り方パターン(2^n通りある)ごとに、その「とり方パターンが真である」という仮説をどれくらいしんじたらよいか、その仮説の場合に特定のphenotypeでの成否確率分布を値ベクトルpについて計算する
# 説明変数の採用パターンごとの尤度の変化の比(すべての説明変数で2^n通りに分けた場合を基準とする)
K <- rep(0,length(X[,1]))
# 説明変数の採用パターンごとの、特定phenotypeでの成否確率分布
prs <- matrix(0,length(K),length(p))
for(i in 1:length(X[,1])){
	# この仮説で問題にする説明変数を取り出す
	s <- which(X[i,]==1)
	#print(s)
	# apply()関数で部分分割表化する
	tmp1 <- my.marginal.apply(A1,s,sum)
	tmp2 <- my.marginal.apply(A2,s,sum)
	# 成否ごとのカウントとその和とが必要
	marg <- tmp1+tmp2
	#print(tmp1)
	#print(tmp2)
	# 限定した仮説に関して、フルモデルのもとでの尤度にもとづく「比」を計算する
	K[i] <- sum(lgamma(marg+2)-lgamma(tmp1+1)-lgamma(tmp2+1))
	# 仮説においては、特定のセルにおける成否カウントとのその和に基づくベータ分布が必要になる
	# apply()関数で作った部分分割表のどのセルがその特定のセルかを与える
	loc <- 1
	if(length(s)>0){
		loc <- sum(phenotype[s]*2^(0:(length(s)-1)))+1
	}
	#print(loc)
	#prs[i,] <- log(dbeta(p,tmp1[length(tmp1)]+1,tmp2[length(tmp2)]+1))
	
	prs[i,] <- dbeta(p,tmp1[loc]+1,tmp2[loc]+1)
	#for(j in 1:length(tmp1)){
	#	prs[i,] <- prs[i,] + log(dbeta(p,tmp1[j]+1,tmp2[j]+1))
	#}
	#print(which(X[i,]==1))
	#print(tmp1)
	#print(tmp1)
}
# フルモデルを基準に比を取る
# 対数で計算していたから、加減計算して、指数関数にかける
K.2 <- -K+K[length(K)]
exp(K.2)
# 仮説の事前確率はすべて同一だったとする
pre.p <- rep(1/length(X[,1]),length(X[,1]))
# これが仮説別の事後確率
post.p <- pre.p * exp(K.2)
post.p <- post.p/sum(post.p)
post.p
# 仮説別の事後確率で重みづけをして、統合的成功事後確率分布を求める
PR <- rep(0,length(p))
for(i in 1:length(pre.p)){
	PR <- PR + post.p[i]*prs[i,]
}
#exp.prs <- exp(prs)

par(mfcol=c(1,2))
matplot(p,cbind(PR,t(exp.prs)),type="l",ylim=c(0,max(exp.prs)),lwd=c(5,rep(1,length(exp.prs))))
plot(p,PR,type="l",ylim=c(0,max(exp.prs)))
par(mfcol=c(1,1))

#matplot(p,cbind(est.women,PR),type="l")