# The data for Thisted's widow example.    #
#   y0=3062 widows had 0 children...       # 
#   y[6]=2  had 6. #

y<-c(3062,587,284,103,33,4,2)


# Putting in the negative log-likelihood   #
#   for the Binomial-Poisson mixture and   #
#   using optim directly.                  #

cnloglike<-function(pars,dat){
 n<-length(dat)
 p<-pars[1]
 m<-pars[2]
 y0<-dat[1]
 y<-dat[2:n]
 cnll<-y0*log(p+(1-p)*exp(-m))
 for (i in 1:(n-1))
   {cnll<-cnll+y[i]*(log(1-p)+i*log(m)-m-sum(log(1:i)))}
 -cnll}

optim(c(.5,1),cnloglike,dat=y) 


# The EM formulation where y0 is partititioned #
#   into the xA who were binomial and the XB   #
#   who were Poisson.                          #  

emwidow<-function(init,dat,niter=20){
 n<-length(dat)
 p<-init[1]
 m<-init[2]
 y0<-dat[1]
 y<-dat[2:n]
 xA<-0
 xB<-0
 for (i in 1:niter){
   xA<-c(xA,y0*p[i]/(p[i]+(1-p[i])*exp(-m[i])))
   xB<-c(xB,y0-xA[i+1])
   p<-c(p,xA[i+1]/(y0+sum(y)))
   m<-c(m,sum((1:(n-1))*y)/(xB[i+1]+sum(y)))
 }
 cbind(xA,xB,p,m)
} 
 
emwidow(c(.5,1),y,niter=60)