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