---
title: "em for recessive allele"
author: "Jerome Goudet"
date: '`r Sys.Date()`'
output: html_document
---

The following code will estimate the recessive allele frequency from the counts of the dominant phenotype assuming Hardy Weinberg (see slides 80-81). 

```{r}
n<-1000
naa<-23
dat<-c(n-naa,naa)
n<-sum(dat)
(pa1<-(dat[2]/n)^.5)

em.ra<-function(dat,ini=NULL,print.it=TRUE,eps=0.0001){
  if (is.null(ini)) pa2<-runif(1) else pa2<-ini 
  dif<-1
  acc<-0
  while (dif>=eps){
    acc<-acc+1
    pa2n<-2*(n*pa2+dat[2])/(2*n*(1+pa2))
    dif<-abs(pa2n-pa2)
    if (print.it) cat(paste(acc," : ",round(pa2,digits=3)," -> ",round(pa2n,digits=3),"\n"))
    pa2<-pa2n
  }
  pa2
}

em.ra(dat)

```