#  R program to calculate maximum likelihood (ML) estimates for parameters
#  a b and o (allele frequencies) in the multinomial model for 
#  Hardy-Weinberg proportions in human blood types.  The model is
#  Y1,Y2,Y3,Y4 ~ multinomial(n,p1,p2,p3,p4) where
#  p1=a+a + 2*a*o  blood type A proportion
#  p2=b*b + 2*b*o  blood type B proportion
#  p3=2*a*b        blood type AB proportion
#  p4=o*o          blood type O proportion
#  
#  Here 0<a<1, 0<b<1, o=1-a-b.

#  Count frequencies are entered into the vector yy here.
yy=c(182,60,17,176)

#  Data in example are from:  Rao, CR 1973. Linear statistical inference and
#  its applications.  Wiley.

#  Set initial patameter values here.
a0=.33
b0=.33

#  ML objective function "negloglike.ml" is negative of log-likelihood;
#  the Nelder-Mead optimization routine in R, "optim", is a minimization
#  routine.  The two function arguments are:  theta = vector of
#  parameters (transformed to real line), ys = vector of frequencies.

negloglike.ml=function(theta,ys)  
{
   a=exp(-exp(theta[1]))     #  Constrains 0 < a < 1.
   b=exp(-exp(theta[2]))     #  Constrains 0 < b < 1.
   o=1-a-b
   n=sum(ys)
   k=length(ys)
   p=c( a*a+2*a*o,
        b*b+2*b*o,
        2*a*b,
        o*o )              #  H-W model for the probabilities.

   ofn=-sum(ys*log(p))     #  No need to calculate all the factorials.
   return(ofn)
}

# The ML estimates.
MULTML=optim(par=c(log(-log(a0)),log(-log(b0))),
   negloglike.ml,NULL,method="Nelder-Mead",ys=yy)
results=c(exp(-exp(MULTML$par[1])),exp(-exp(MULTML$par[2])),-MULTML$val)
a.ml=results[1]            # These are the ML estimates.
b.ml=results[2]            #          --
o.ml=1-a.ml-b.ml           #          --

nn=sum(yy)
loglike.ml=results[3]+lfactorial(nn)-sum(lfactorial(yy))

# Calculate expected values, LR statistic, etc.
pp=c( a.ml*a.ml+2*a.ml*o.ml,
        b.ml*b.ml+2*b.ml*o.ml,
        2*a.ml*b.ml,
        o.ml*o.ml )
EE=nn*pp
y1=yy
y1[y1==0]=1                # Guard against log(0) in G-squared.
Gsq=2*sum(yy*log(y1/EE))   # G-squared goodness of fit statistic.
pval=1-pchisq(Gsq,1)
Xsq=sum((yy-EE)^2/EE)

#  Print the results.
a.ml
b.ml
o.ml
loglike.ml
Gsq
pval
Xsq
cbind(EE,yy)