EM - algorithm: LD

### LD
library(genetics)
fms<-read.delim(file="http://people.umass.edu/foulkes/asg/data/FMS_data.txt",header=T,sep="\t")
attach(fms)

Actn3Snp1<-genotype(actn3_r577x,sep="")
Actn3Snp2<-genotype(actn3_rs540874,sep="")
Actn3Snp3 <- genotype(actn3_rs1815739,sep="")
Actn3Snp4 <- genotype(actn3_1671064,sep="")

LD(Actn3Snp2,Actn3Snp3)

tt<-table(Actn3Snp2,Actn3Snp3)
temp<-tt[1,]
tt[1,]<-tt[3,];tt[3,]<-temp

pAB<-0.25;paB<-0.25;pAb<-0.25;pab<-0.25

N<-sum(tt)
pAB<-(2*tt[1,1]+tt[1,2]+tt[2,1]+tt[2,2]*(pAB*pab)/(pAB*pab+paB*pAb))/N/2
pAb<-(2*tt[1,3]+tt[1,2]+tt[2,3]+tt[2,2]*(paB*pAb)/(pAB*pab+paB*pAb))/N/2
paB<-(2*tt[3,1]+tt[3,2]+tt[2,1]+tt[2,2]*(paB*pAb)/(pAB*pab+paB*pAb))/N/2
pab<-(2*tt[3,3]+tt[3,2]+tt[2,3]+tt[2,2]*(pAB*pab)/(pAB*pab+paB*pAb))/N/2
cbind(pAB,pAb,paB,pab); pAB+pAb+paB+pab; 
2*tt[1,1]+tt[1,2]+tt[2,1]+tt[2,2]*(pAB*pab)/(pAB*pab+paB*pAb)+
2*tt[1,3]+tt[1,2]+tt[2,3]+tt[2,2]*(paB*pAb)/(pAB*pab+paB*pAb)+
2*tt[3,1]+tt[3,2]+tt[2,1]+tt[2,2]*(paB*pAb)/(pAB*pab+paB*pAb)+
2*tt[3,3]+tt[3,2]+tt[2,3]+tt[2,2]*(pAB*pab)/(pAB*pab+paB*pAb)
tt[1,1]<-pAB^2
tt[1,2]<-2*pAB*pAb
tt[1,3]<-pAb^2
tt[2,1]<-2*pAB*paB
tt[2,2]<-2*pAB*pab+2*paB*pAb
tt[2,3]<-2*pAb*pab
tt[3,1]<-paB^2
tt[3,2]<-2*pab*paB
tt[3,3]<-pab^2
tt<-tt*N;tt

row<-apply(tt,1,sum); col<-apply(tt,2,sum)

pA<-(2*row[1]+row[2])/sum(row)/2;
pB<-(2*col[1]+col[2])/sum(col)/2;

pAB-pA*pB

EM - algorithm: ABO

### ABO
Na<-25; Nab<-50; Nb<-25; No<-15; N<-Na+Nb+No+Nab

pa<-0.33; pb<-0.33; po<-0.34

for(i in 1:20){
Naa<-pa^2/(pa^2+2*pa*po)*Na
Nao<-2*pa*po/(pa^2+2*pa*po)*Na
Nbb<-pb^2/(pb^2+2*pb*po)*Nb
Nbo<-2*pb*po/(pb^2+2*pb*po)*Nb
#Nab<-2*pa*pb*Nab
#Noo<-po^2*No
pa<-(2*Naa+Nab+Nao)/(Naa+Nao+Nbb+Nbo+Nab+No)/2;
pb<-(2*Nbb+Nab+Nbo)/(Naa+Nao+Nbb+Nbo+Nab+No)/2;
po<-(2*No+Nao+Nbo)/(Naa+Nao+Nbb+Nbo+Nab+No)/2;
}

cbind(Naa,Nao,Nbb,Nbo,Nab,No)
cbind(pa,pb,po)