#The vector of counts, y, should be ordered so that the first variable #has the fastest changing index, with 1 appearing first, then 0. #The sufficient statistics are inefficient to compute in every iteration, #so they're included as arguments. See the code at the end of #margcon to compute sufficient statistics and an initial estimate for beta margcon=function(betahat,y,indx,ypt,scor,ns){ fullbeta=c(betahat,0) #Numerator terms num=sum(ypt*fullbeta) #Denominator terms expb=exp(indx%*%fullbeta) sumscor=tapply(expb,scor,sum) den=sum(ns*log(sumscor)) logL=den-num return(logL) } #Data y=c(571,11,64,19,67,18,18,30,66,14,36,39,63,53,49,270) ny=length(y) T=log2(ny) #Create table indices to compute sufficient statistics for beta indx=NULL for(i in 1:T){ indx=cbind(indx,(rep(rep(c(1,0),each=2^(i-1)),2^(T-i)))) } ypt=apply(indx*y,2,sum) scor=apply(indx,1,sum) ns=tapply(y,scor,sum) #Initial values for beta #Scale estimates between -2 and 2 betahat=2*(ypt-ypt[T])/(max(abs(ypt-ypt[T]))) #beta[T] is fixed at 0 betahat=betahat[-T] nlminb(betahat,margcon,y=y,indx=indx,ypt=ypt,scor=scor,ns=ns) #Use optim to compute standard errors betaout=optim(betahat,margcon,method="BFGS",hessian=T,y=y,indx=indx,ypt=ypt,scor=scor,ns=ns) sqrt(diag(solve(betaout$hessian))) #Compute ratio statistics expb=exp(indx%*%c(betaout$par,0)) sumscor=tapply(expb,scor,sum) #Confirm that the ratio statistics below form a moment sequence for a non-negative measure ratiostat=ns/sumscor mommatrix=matrix(c(ratiostat[1:3],ratiostat[2:4],ratiostat[3:5]),ncol=3) #All the eigenvalues of the following two symmetric matrices are positive, #hence both matrices are positive definite, and our conditional MLEs are #the same as the joint MLEs for the random effects model eigen(mommatrix[1:3,1:3]) eigen(mommatrix[2:3,1:2])