DNA鑑定と幾何平均

  • こちらでDNA鑑定とバースデイ・パラドクスの話を書いた
  • その延長線上で、同時確率の「期待値」を考えたいのだが、「期待値」が算術平均なのに対して、幾何平均なら簡単に計算できるのだが、幾何平均を持って、「代表値」とすることの意義はあるのか、あるとしたらどのあたりにあるのか…はまだ不明
# AmpFlSTRIdentifilerの15 STR

Identifiler15<-c("D8S1179","D21S11","D7S820","CSF1PO","D3S1358","TH01","D13S317","D16S539","D2S1338","D19S433","vWA","TPROX","D18S51","D5S818","FGA")

AllelesId<-NULL
ProbsIDYoshidaN<-NULL
ProbsIDYoshidaFreq<-NULL
ProbsIDYoshidaFreq2<-NULL
AllelesId[[1]]<-c("7",9,10,11,12,13,14,15,16,17,18)
AllelesId[[2]]<-c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2)
AllelesId[[3]]<-c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15)
AllelesId[[4]]<-c("7",8,9,10,11,12,13,14,15,16)
AllelesId[[5]]<-c("12",13,14,15,16,17,18,19,21)
AllelesId[[6]]<-c("5",6,7,8,9,9.3,10)
AllelesId[[7]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[8]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[9]]<-c("16",17,18,19,20,21,22,23,24,25,26,27,28)
AllelesId[[10]]<-c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18)
AllelesId[[11]]<-c("13",14,15,16,17,18,19,20,21,22)
AllelesId[[12]]<-c("7",8,9,10,11,12,13,14)
AllelesId[[13]]<-c("10",11,12,13,14,15,16,17,17.1,18,19,20,21,22,23,24,25,26,27)
AllelesId[[14]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[15]]<-c("17",18,19,20,21,22,22.2,23,23.2,24,24.2,25,25.2,26,27,28)

ProbsIDYoshidaN[[1]]<-c(1,2,358,294,331,607,553,363,173,17,1)
ProbsIDYoshidaN[[2]]<-c(3,114,13,666,910,12,3,279,192,52,321,7,6,109,2,11)
ProbsIDYoshidaN[[3]]<-c(8,342,123,1,591,1,2,887,634,94,16,1)
ProbsIDYoshidaN[[4]]<-c(27,3,126,603,561,1130,186,48,14,2)
ProbsIDYoshidaN[[5]]<-c(5,2,78,1069,827,539,171,8,1)
ProbsIDYoshidaN[[6]]<-c(4,603,720,179,1076,94,24)
ProbsIDYoshidaN[[7]]<-c(3,721,348,310,597,546,138,35,2)
ProbsIDYoshidaN[[8]]<-c(1,5,955,531,505,482,196,22,3)
ProbsIDYoshidaN[[9]]<-c(23,263,429,564,285,40,136,396,291,166,78,23,6)
ProbsIDYoshidaN[[10]]<-c(1,3,10,2,109,14,776,81,944,238,137,310,14,52,8,1)
ProbsIDYoshidaN[[11]]<-c(1,524,72,497,763,609,199,27,7,1)
ProbsIDYoshidaN[[12]]<-c(1,1227,308,83,980,96,3,2)
ProbsIDYoshidaN[[13]]<-c(6,12,129,538,600,454,339,220,1,130,98,59,42,35,20,10,4,2,1)
ProbsIDYoshidaN[[14]]<-c(7,18,232,542,789,635,450,24,3)
ProbsIDYoshidaN[[15]]<-c(9,58,181,240,353,544,5,554,13,424,3,198,5,86,21,6)
for(i in 1:15){
	ProbsIDYoshidaFreq[[i]]<-ProbsIDYoshidaN[[i]]/sum(ProbsIDYoshidaN[[i]])
	ProbsIDYoshidaFreq2[[i]]<-(ProbsIDYoshidaN[[i]]+1)/(sum(ProbsIDYoshidaN[[i]])+1)
	ProbsIDYoshidaFreq2[[i]][which(ProbsIDYoshidaN[[i]]<5)]<-5/sum(ProbsIDYoshidaN[[i]])
}



# マーカー座位数
n.m <- 15
# マーカー別アレル数
#n.a <- sample(4:20,n.m,replace=TRUE)
# マーカー別アレル頻度(対数)
library(MCMCpack)
A <- list()
for(i in 1:n.m){
	#A[[i]] <- c(rdirichlet(1,rep(1,n.a[i])))
	A[[i]] <- ProbsIDYoshidaFreq2[[i]]
	A[[i]] <- A[[i]]/sum(A[[i]])
}

# マーカー別ジェノタイプ頻度と
# マーカー別ジェノタイプ種類数
D <- list()
n.d <- rep(0,n.m)
for(i in 1:n.m){
	tmp <- outer(A[[i]], A[[i]], FUN="*")
	tmp <- tmp + t(tmp)
	diag(tmp) <- diag(tmp)/2
	D[[i]] <- log(c(diag(tmp),tmp[which(upper.tri(tmp))]))
	
	n.d[i] <- length(D[[i]])
}

A
D

lapply(A,sum)
lapply(lapply(D,exp),sum)

# 複合ジェノタイプ数
L.whole.log <- sum(log(n.d))

# マーカー別に、当該マーカー以外のマーカーが作るジェノタイプ数
L.marker.log <- log(n.d)

L.marker.log

# 人数
#k <- sample((exp(L.whole.log)/2):exp(L.whole.log),1)

k <- log(1.25 * 10^8)
# 複合ジェノタイプから重複なしにk人分のジェノタイプを選んだ時に
# そのジェノタイプの確率の積がいくつになるかを
# 幾何平均として求める

# 最頻・最稀複合ジェノタイプ
max.g <- 0
min.g <- 0
for(i in 1:n.m){
	max.g <- max.g + max(D[[i]])
	min.g <- min.g + min(D[[i]])
}
max.g
min.g

exp(max.g)
exp(min.g)

ret <- 0
geom.g <- 0

for(i in 1:n.m){
	ret <- ret + sum(D[[i]]) * exp(-L.marker.log[i]+L.whole.log+k-L.whole.log)
	geom.g <- geom.g + sum(D[[i]]) * exp(-L.marker.log[i]+L.whole.log-L.whole.log)
	#geom.g <- geom.g + sum(D[[i]]) * exp(-L.marker.log[i])
	#ret <- ret + sum(D[[i]]) * exp(-L.marker.log[i] + L.whole.log) * k/exp(L.whole.log)
}

geom.g
ret

exp(geom.g)

exp(ret)