selfing sucks

#make a random haplotype given a vector of allele freqs.
ran.hap=function(numloci,p){sapply(1:numloci,function(x) rbinom(1,1,p[x]))}

#make a kid from two parents. simple version with no recombination. includes het and hom error
kid=function(mom,dad,het.error,hom.error){
    k1=copy.mom(mom)
    k2=copy.mom(dad)

    truekid=k1[[2]]+k2[[2]]
    hets_with_error=sample(which(truekid==1),round(het.error*length(which(truekid==1))))
    hom0_with_error=sample(which(truekid==0),round(hom.error*length(which(truekid==0))))
    hom2_with_error=sample(which(truekid==2),round(hom.error*length(which(truekid==2))))
    truekid=replace(truekid,hets_with_error,sample(c(0,2),length(hets_with_error),replace=T)  )
    truekid=replace(truekid,hom0_with_error,sample(c(1,2),length(hom0_with_error),replace=T)  )
    truekid=replace(truekid,hom2_with_error,sample(c(1,0),length(hom2_with_error),replace=T)  )
    return(list(truekid,k1[[1]]+k2[[1]]))
}

#ML estimation of three genotypes
probs<-function(e,eprime,j,m,n){
  phet=((m+n)*log(e/2)+j*log(1-e))*2; #p of hets. times two is because Mendel 1:2:1
  p00=(n+j)*log(eprime/2)+m*log(1-eprime); #p of 00 hom
  p11=(m+j)*log(eprime/2)+n*log(1-eprime); #p of 11 hom

  return(c(p00,phet,p11)) # positive numbers mean het more likely
}

#simple version of copying hapltoypes from mom, with no recombination
copy.mom=function(mom){  
    momhap=rbinom(1,1,.5)+1
    k1=mom[[momhap]]    
    return(list(momhap,k1))
}

params

numloci=1000
het.error=0.7
hom.error=0.01

#set up SFS
x=1:99/100
freq=1/x
sfs=as.numeric(); 
for(i in 1:99){sfs=c(sfs,rep(x[i],100*freq[i]))}

make mom,selfed kid

  wins=0
for(i in 1:100){
  p=sample(sfs,numloci)
  a1=ran.hap(numloci,p)
  a2=ran.hap(numloci,p)
  genos<-which(a1!=a2)
  mom=list(rep(0,length(genos)),rep(1,length(genos)))
  #dude=kid(mom,mom,p,0,min(het.error*(rnorm(1,1,0.2)),1),hom.error*(rnorm(1,1,0.2)))
  dude=kid(mom,mom,0.7,0.002)
  hets=length(which(dude[[1]]==1))
  hom0=length(which(dude[[1]]==0))
  hom1=length(which(dude[[1]]==2))
  logprobs=probs(het.error,hom.error,hets,hom0,hom1)
  if(logprobs[dude[[2]]-1]==max(logprobs)){ wins=wins+1;}
#  if(i %% 100 ==0){ print(i)}
}
print(wins)

## [1] 100