一致確率

• マーカー数
• 第iマーカーのアレル数
• 各マーカーのアレル頻度
• 各マーカーのディプロタイプ数
• 各マーカーのディプロタイプ頻度
  • HWEを仮定すれば、はから導ける
• あるマーカーでディプロタイプが一致する確率は、
• すべてのマーカーでディプロタイプが一致する確率は

# allele freq ベクトルからdiplotype freqベクトルを作る
diplotypeFreq<-function(fa){
	m<-outer(fa,fa,"*")
	m2<-2*m
	diag(m2)<-diag(m2)/2
	m2[!lower.tri(m2)]
}
# diplotype freqベクトルから、全マーカー一致確率の自然対数値を計算する
ProbIdenticalLN<-function(ds){
	tmp<-0
	for(i in 1:length(ds)){
		p<-sum(ds[[i]]^2)
		if(p==0){
			return(-Inf)
		}else{
			tmp<-tmp+log(p)
		}
	}
	return(tmp)
}

# HWE下
# アレル頻度given
# diplotype 一致確率
# diplotype freqを作ろう
Nm<-15
library(MCMCpack)
Nas<-sample(2:10,Nm,replace=TRUE)
fs<-list()
ds<-list()
for(i in 1:Nm){
	fs[[i]]<-c(rdirichlet(1,rep(1,Nas[i])))
	print(diplotypeFreq(fs[[i]]))
	ds[[i]]<-diplotypeFreq(fs[[i]])
}
# 一致確率を計算しよう
ProbIdenticalLN(ds)