メイティングして生じるディプロタイプ頻度を計算する２

昨日の記事の続き
phase knownな記載にしよう
昨日の記事はphase unknownで、上三角行列(と対角成分)にしたところを、正方行列全体を使う
アレル数 $N_a$ の多型では、個人のディプロタイプを母由来を行、父由来を列とすれば、 $N_a\times N_a$ の正方行列で表される

Na<-5
MakeRandomDiplotypeFreq<-function(Na,af){
 p<-runif(Na^2)
 p<-p/sum(p)
 matrix(p,Na,Na)
}
P<-MakeRandomDiplotypeFreq(Na)
P
sum(P)

この個人の配偶子のアレル頻度を求めてみよう

gameteFreq<-function(P){
 (apply(P,1,sum)+apply(P,2,sum))/2
}
G<-gameteFreq(P)
G
sum(G)
||
-こんな個人２人から生まれる子は
>|r|
OffspringDiplotypeFreq<-function(m,f){# m:mother,f:father
 gameteM<-gameteFreq(m)
 gameteF<-gameteFreq(f)
 outer(gameteM,gameteF,FUN="*")
}
m<-MakeRandomDiplotypeFreq(Na)
f<-MakeRandomDiplotypeFreq(Na)
o<-OffspringDiplotypeFreq(m,f)
o
sum(o)

父母のディプロタイプが決まっていれば

RandomDiplotypeMatrix<-function(Na){
 m<-matrix(0,Na,Na)
 m[sample(1:(Na^2),1)]<-1
 m
}
m<-RandomDiplotypeMatrix(Na)
f<-RandomDiplotypeMatrix(Na)
m
f
o<-OffspringDiplotypeFreq(m,f)
o

父　母　子のディプロタイプがあれば、尤度は。。。
- 父・母ともにホモのとき、子のディプロタイプは１パターンしか出ない。
- 父・母・子のディプロタイプを勝手に振って、そのようになる場合がなくはないことが以下からわかる

Na<-3
Niter<-10000
res<-rep(0,Niter)
for(i in 1:Niter){
m<-RandomDiplotypeMatrix(Na)
f<-RandomDiplotypeMatrix(Na)
o<-RandomDiplotypeMatrix(Na)
res[i]<-sum(OffspringDiplotypeFreq(m,f)*o)
}
hist(res)

じゃあ、父、母のディプロタイプが(確率的に)与えられたときに、その子のディプロタイプ頻度の期待値が行列演算で出るから、その期待頻度からの乱数を発生させることで、子供をたくさん産んでみよう

Na<-4
m<-RandomDiplotypeMatrix(Na)
f<-RandomDiplotypeMatrix(Na)

Nkids<-100

kidsFreq<-OffspringDiplotypeFreq(m,f)

Sibs<-sample(1:(Na^2),Nkids,prob=c(kidsFreq),replace=TRUE)
SibsVector<-rep(0,Na^2)
for(i in 1:length(Sibs)){
	SibsVector[Sibs[i]]<-SibsVector[Sibs[i]]+1
}
SibsMatrix<-matrix(SibsVector,Na,Na)

SibsFreq<-SibsMatrix/sum(SibsMatrix)
SibsFreq
kidsFreq