normIntFunc=function (u,b,sc,sig) {
	tb=length(b)
	odds=exp(sig*u*sc)
	for (i in 1:tb) {
		odds=odds/(1+exp(sig*u+b[i]))
	}
	sam = odds * dnorm(u)
	return(sam)
}

logL=function (parm,scor,indx,y) {
S=length(parm)
T=S-2
ny=length(y)
beta=parm[1:T]
sig=parm[S-1]
ntot=parm[S]
piI0=NULL
for(i in 1:ny){
w1=integrate(normIntFunc,-4,4,b=beta,sc=scor[i],sig=sig)$value
piI0[i]=w1*as.numeric(exp(indx[i,]%*%beta))
}
piI01=integrate(normIntFunc,-4,4,b=beta,sc=0,sig=sig)$value
n0=ntot-sum(y)
lf=as.numeric(lfactorial(ntot)-sum(lfactorial(y))-lfactorial(n0)+y %*% log(piI0)+n0*log(piI01))
return(-lf)
}

#Compute statistics
y=c(3,6,0,5,1,0,0,3,2,3,0,0,1,0,0,4,2,3,1,0,1,0,0,1,0,0,0,0,0,0,0,4,1,1,1,2,0,2,0,4,0,3,0,1,0,2,0,2,0,1,0,1,0,1,0,1,1,1,0,0,0,1,2)
ny=length(y)
T=log2(ny+1)
indx=NULL
#Construct response vector
for(i in 1:T){
indx=cbind(indx,rep(rep(c(0,1),each=2^(i-1)),2^(T-i)))
}
indx=indx[-1,]
sigmaHat = 1.0
betaHat=2*(ymarg-ymarg[T])/max(abs(ymarg-ymarg[T]))
ntot=90
parmHat = c(betaHat, sigmaHat,ntot)
#Column marginals
ymarg=apply(indx*y,2,sum)
#Compute score for each cell
scor=apply(indx,1,sum)
#Score group totals
ns=tapply(y,scor,sum)
#Check loglikelihood function
logL(parmHat,scor,indx,y)
#Maximize
soln=nlminb(start=parmHat,logL,lower=c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,0,0),scor=scor,indx=indx,y=y)

#CI
LogL=NULL
for(ntot in 75:110){
parmHat=c(soln$par[1:T],soln$par[T+1],ntot)
LogL[ntot-74]=2*logL(parmHat,scor,indx,y)
}
refline=2*logL(soln$par,scor,indx,y)+3.84
plot(75:110,LogL,type="l")
abline(h=refline,col="red")
#We could use a root-finder to obtain endpoints, but the arguments for logL would have to #be rewritten