- SNPのディプロタイプ頻度をとして表したときに、観察3ディプロタイプ人数を用いてその対数尤度はと表される
- HWEを仮定すると変数fは0に固定することになるので、変数はpのみ
- HWEの仮定とHWDの仮定では変数の個数が2,1と、その差が1。このためHWE検定の自由度は1
- 今適当に、p,fを与えて、Gを作り、LLを最大にするp,fをグリッド計算で出してみよう
p<-0.3
f<-0.1
N<-100
g<-c(p^2+p*(1-p)*f,2*p*(1-p)*(1-f),(1-p)^2+p*(1-p)*f)*N
Ps<-seq(from=0.01,to=0.99,by=0.01)
fs<-seq(from=-0.99,to=0.99,by=0.01)
Pfs<-expand.grid(Ps,fs)
Ls2<-rep(-Inf,length(Pfs[,1]))
for(i in 1:length(Pfs[,1])){
tmpp<-Pfs[i,1]
tmpf<-Pfs[i,2]
tmpG<-c(tmpp^2+tmpp*(1-tmpp)*tmpf,2*tmpp*(1-tmpp)*(1-tmpf),(1-tmpp)^2+tmpp*(1-tmpp)*tmpf)
if(min(tmpG>0)){
Ls2[i]<-sum(g*log(tmpG))
}
}
library(rgl)
Ls2[which(Ls2==-Inf)]<-min(Ls2[which(Ls2!=-Inf)])
plot3d(Pfs[,1],Pfs[,2],Ls2)
maxpoint<-which(Ls2==max(Ls2))
Pfs[maxpoint,]