2013-09-01
関数にしてみる:法廷用ベイジアンネットワーク

library(gRain)

my.forensic.bn.2 <- function(n.p.type,N.p,n.g.type,Freq.g,frac.s,typing.prec.mat,L.prob,mot.levels,M.prob,Cr.prob,rs){
	# タイプ別個人の存在
	ps <- Gs <- list()
	yn <- c("1","0")
	p.type <- paste("p",1:n.p.type,sep="")
	for(i in 1:n.p.type){
		ps[[i]] <- cptable(~p,values=c(N.p[i],sum(N.p)-N.p[i]),levels=yn)
		ps[[i]]$vpa[1] <- paste("p",i,sep="")
	}
	# ジェノタイプの事前確率
	G.type <- paste("Gn",1:n.g.type,sep="")
	for(i in 1:n.p.type){
		Gs[[i]] <- cptable(~G,values=Freq.g,levels=G.type)
		Gs[[i]]$vpa[1] <- paste("G",i,sep="")
	}
	# 容疑者のサンプルをタイピング
	# 取り違えずに標本を取る
	tmp.exp <- expand.grid(rep(list(1:n.g.type),n.p.type))
	tmp.val <- matrix(0,length(tmp.exp[,1]),n.g.type)
	for(i in 1:length(tmp.val[,1])){
		tmp.val[i,tmp.exp[i,1]] <- tmp.val[i,tmp.exp[i,1]]+frac.s
		for(j in 2:n.p.type){
			tmp.val[i,tmp.exp[i,j]] <- tmp.val[i,tmp.exp[i,j]]+(1-frac.s)/(n.p.type-1)
		}
	}
	Sx <- cptable(~Sx,values=c(t(tmp.val)),levels=G.type)
	for(i in 1:n.p.type){
		Sx$vpa <- c(Sx$vpa,Gs[[i]]$vpa[1])
	}
	#tmp.mat <- matrix((1-precision.typing)/(n.g.type-1),n.g.type,n.g.type)
	#diag(tmp.mat) <- precision.typing
	# タイピングエラーを考慮
	SGx <- cptable(~SGx|Sx,values=c(t(typing.prec.mat)),levels=G.type)

	# 容疑者標本のbn
	tmplist <- list(Sx,SGx)
	tmplist1 <- list(Sx,SGx)
	for(i in 1:length(ps)){
		tmp.n <- length(tmplist)
		tmplist[[tmp.n+1]] <- ps[[i]]
		tmplist[[tmp.n+2]] <- Gs[[i]]
		tmplist1[[length(tmplist1)+1]] <- Gs[[i]]
	}
	#plist1 <- compileCPT(tmplist1)
	#net1 <- grain(plist1)

	# 現場への出入り
	
	Ls <- list()
	for(i in 1:n.p.type){
		Ls[[i]] <- cptable(~l,values=c(L.prob[i,],0,1),levels=yn)
		Ls[[i]]$vpa[1] <- paste("L",i,sep="")
		Ls[[i]]$vpa <- c(Ls[[i]]$vpa,ps[[i]]$vpa[1])
		tmplist[[length(tmplist)+1]] <- Ls[[i]]
	}
	#plist2 <- compileCPT(tmplist)
	#net2 <- grain(plist2)

	# 動機
	Ms <- list()
	for(i in 1:n.p.type){
		Ms[[i]] <- cptable(~m,values=M.prob[i,],levels=mot.levels)
		Ms[[i]]$vpa[1] <- paste("M",i,sep="")
		Ms[[i]]$vpa <- c(Ms[[i]]$vpa,Ls[[i]]$vpa[1])
		tmplist[[length(tmplist)+1]] <- Ms[[i]]
	}
	#plist3 <- compileCPT(tmplist)
	#net3 <- grain(plist3)

	# 犯行は現場に居ないと犯せない
	# 現場に居て、動機があっても犯すとは限らない
	# 現場に居ると、動機があってもなくても、何かの拍子に犯すかも
	Cs <- list()
	for(i in 1:n.p.type){
		tmp <- cbind(Cr.prob,1-Cr.prob)
		Cs[[i]] <- cptable(~cr,values=c(t(tmp),rep(0:1,length(mot.levels))),levels=yn)
		Cs[[i]]$vpa[1] <- paste("C",i,sep="")
		Cs[[i]]$vpa <- c(Cs[[i]]$vpa,Ms[[i]]$vpa[1],Ls[[i]]$vpa[1])
		tmplist[[length(tmplist)+1]] <- Cs[[i]]
	}
	
	# n.p.type人のうち１人しか犯行は犯せない
	tmp.val <- expand.grid(rep(list(1:0),n.p.type))
	ss <- apply(tmp.val,1,sum)
	ari <- which(ss==1)
	tmp.val2 <- rep(0,length(ss))
	tmp.val2[ari] <- 1
	tmp.val3 <- 1-tmp.val2
	#Cr <- cptable(~Cr,values=c(rep(1:0,2^n.p.type-1),0,1),levels=yn)
	Cr <- cptable(~Cr,values=c(t(cbind(tmp.val2,tmp.val3))),levels=yn)
	for(i in 1:n.p.type){
		Cr$vpa <- c(Cr$vpa,Cs[[i]]$vpa[1])
	}
	tmplist[[length(tmplist)+1]] <- Cr
	#plist4 <- compileCPT(tmplist)
	#net4 <- grain(plist4)

	# 犯人は誰だ
	tmp.mat <- expand.grid(rep(list(1:0),n.p.type))
	sum.tmp.mat <- apply(tmp.mat,1,sum)
	tmp.mat[which(sum.tmp.mat==0),] <- 1
	sum.tmp.mat <- apply(tmp.mat,1,sum)
	tmp.mat2 <- tmp.mat/sum.tmp.mat
	E <- cptable(~E,values=c(t(tmp.mat2)),levels=p.type)
	for(i in 1:n.p.type){
		#E$vpa <- c(E$vpa,Cs[[n.p.type+1-i]]$vpa[1])
		E$vpa <- c(E$vpa,Cs[[i]]$vpa[1])
	}
	tmplist[[length(tmplist)+1]] <- E
	#plist5 <- compileCPT(tmplist)
	#net5 <- grain(plist5)

	# 現場に居たか居ないか、犯行に関与したかしないかで遺留品を残す確率は変わる
	Rs <- list()
	for(i in 1:n.p.type){
		tmp.val <- c(rep(c(rs[2],1-rs[2]),n.p.type),rep(c(0,1),n.p.type))
		tmp.val[c((i-1)*2+1,i*2)] <- c(rs[1],1-rs[1])
		Rs[[i]] <- cptable(~r,values=tmp.val,levels=yn)
		Rs[[i]]$vpa[1] <- paste("R",i,sep="")
		Rs[[i]]$vpa <- c(Rs[[i]]$vpa,E$vpa[1],Ls[[i]]$vpa[1])
		tmplist[[length(tmplist)+1]] <- Rs[[i]]
	}
	#plist6 <- compileCPT(tmplist)
	#net6 <- grain(plist6)
	# n.p.type人のうち１人は「残す」
	Rr <- cptable(~Rr,values=c(rep(1:0,2^n.p.type-1),0,1),levels=yn)
	for(i in 1:n.p.type){
		Rr$vpa <- c(Rr$vpa,Rs[[i]]$vpa[1])
	}
	tmplist[[length(tmplist)+1]] <- Rr

	# 遺留品が誰に帰属するか
	A <- cptable(~A,values=c(t(tmp.mat2)),levels=p.type)
	for(i in 1:n.p.type){
		A$vpa <- c(A$vpa,Rs[[i]]$vpa[1])
	}
	tmplist[[length(tmplist)+1]] <- A
	#plist7 <- compileCPT(tmplist)
	#net7 <- grain(plist7)


	# 遺留品の真のジェノタイプ
	
	exx <- expand.grid(rep(list(1:n.g.type),n.p.type+1))
	exx.1 <- exx[,1:n.p.type]
	exx.2 <- exx[,n.p.type+1]
	tmp.mat <- matrix(0,length(exx.2),n.g.type)
	for(i in 1:length(exx.2)){
		#tmp.mat[i,exx.1[i,n.p.type-exx.2[i]+1]] <- 1
		tmp.mat[i,exx.1[i,exx.2[i]]] <- 1
	}
	Ax <- cptable(~Ax,values=c(t(tmp.mat)),levels=G.type)
	for(i in 1:n.p.type){
		#Ax$vpa <- c(Ax$vpa,Gs[[n.p.type-i+1]]$vpa)
		Ax$vpa <- c(Ax$vpa,Gs[[i]]$vpa)
	}
	Ax$vpa <- c(Ax$vpa,A$vpa[1])
	tmplist[[length(tmplist)+1]] <- Ax
	#plist9 <- compileCPT(tmplist)
	#net9 <- grain(plist9)

	#遺留品のタイピング結果
	#tmp.mat <- matrix((1-precision.typing)/(n.g.type-1),n.g.type,n.g.type)
	#diag(tmp.mat) <- precision.typing

	AGx <- cptable(~AGx|Ax,values=c(t(typing.prec.mat)),levels=G.type)
	tmplist[[length(tmplist)+1]] <- AGx
	#plist10 <- compileCPT(tmplist)
	#net10 <- grain(plist10)
	
	plist <- compileCPT(tmplist)
	plist1 <- compileCPT(tmplist1)
	return(list(full.list=plist,sample.list=plist1))
	#return(list(full.list=plist))

}
# 容疑者のタイプ数
n.p.type <- 2
# タイプ別の人数(それぞれ１人,100万人)
N.p <- c(1,1000000)
# 問題にするジェノタイプは２種類
n.g.type <- 2
# その集団頻度
Freq.g <- c(10^(-10),0)
Freq.g[2] <- 1-Freq.g[1]
# サンプルの取り違えはないとする
frac.s <- 1
# タイピング実験は完璧
typing.prec.mat <- matrix(c(1,0,0,1),byrow=TRUE,2,2)
# 容疑タイプ別の現場存在確率。容疑者はアリバイがあり
# 「雑踏での事件」なので、かなりの人がそこにいた可能性がある
L.prob <- matrix(c(0.1,0.9,10000,990000),byrow=TRUE,ncol=2)
# 動機は２種類、ありかなしか
mot.levels <- c("1","0")
# 容疑者には強い動機があり、一般人には現場で動機が発生する可能性がある
M.prob <- rbind(c(1,0,1,0),c(0.01,0.99,0,1))

# 動機と犯行の関係は、「動機があれば、『必ず』犯行に及ぶ」という想定
Cr.prob <- c(1,0)

# 犯行時と非犯行時の遺留率。現場に居れば、必ず遺留品を残し、それは犯行に関わったか否かによらない、とする
rs <- c(1,1)
yn <- c("1","0")

# カテゴリレベルの表記
yn <- c("1","0")
p.type <- paste("p",1:n.p.type,sep="")
G.type <- paste("Gn",1:n.g.type,sep="")

# ベイジアンネットワークを作る
my.bn <- my.forensic.bn.2(n.p.type,N.p,n.g.type,Freq.g,frac.s,typing.prec.mat,L.prob,mot.levels,M.prob,Cr.prob,rs)
my.net.full <- grain(my.bn$full.list)

querygrain(my.net.full,nodes=c("Cr"))
querygrain(my.net.full,nodes=c("C1","C2"),type="joint")
querygrain(my.net.full,nodes=c("E"))