#  R program to calculate profile log-likelihood for population size 
#  using data from a two-sample capture recapture survey.

#  Counts are entered into the vector yy here:  capture histories (1 1),
#  (1 0), (0 1), in that order.
yy=c(7,80,7)

#  Data in example are from: Skalski et al. 1983 Ecology 64:752-760. 

#  Initial patameter values are calculated from multinomial ml estimates.
theta10=yy[1]/(yy[1]+yy[3])
theta20=yy[1]/(yy[1]+yy[2])
psi10=log(theta10/(1-theta10))
psi20=log(theta20/(1-theta20))
nn=sum(yy)
tt0=nn/(1-(1-theta10)*(1-theta20))

#  Range of values for profile plot.  Adjust "tlo" and "thi" to get a
#  good plot.  
tlo=110
thi=370
t.vals=tlo:thi
loglike.vals=t.vals

#  ML objective function "negloglike.ml" is negative of log-likelihood;
#  the optimization routine in R, "optim", is a minimization
#  routine.  The three function arguments are:  psi = vector of
#  parameters (transformed to real line), ts = fixed population size for
#  profile calculation, ys = vector of frequencies.

negloglike.ml=function(psi,ts,ys)  {
   theta1=1/(1+exp(-psi[1]))     #  Constrains 0 < theta1 < 1.
   theta2=1/(1+exp(-psi[2]))     #  Constrains 0 < theta2 < 1.
   k=length(ys)+1
   p=rep(0,k)
   p[1]=theta1*theta2
   p[2]=theta1*(1-theta2)
   p[3]=(1-theta1)*theta2
   p[4]=1-(p[1]+p[2]+p[3])
   nn=sum(ys)
   ys=c(ys,ts-nn)
   ofn=-(lfactorial(ts)-sum(lfactorial(ys))+sum(ys*log(p)))
   return(ofn)
}

# Maximize the log-likelihood function for every value of population size
# from tlo to thi, using the optim() function.  Store each log-likelihood
# value in the vector called loglike.vals.

for (tt in tlo:thi) {
   ii=tt-tlo+1
   MULTML=optim(par=c(psi10,psi20),
      negloglike.ml,NULL,method="Nelder-Mead",ts=tt,ys=yy)
   loglike.vals[ii]=-MULTML$val
}

# Plot profile likelihood and chisquare-based 95% reference line.
ci.height=max(loglike.vals)-3.843/2  #  3.843 = 95th chisquare(1) percentile
ci.cut=rep(ci.height,length(t.vals))
plot(t.vals,loglike.vals,type="l",xlab="population size",
   ylab="profile log-likelihood")
points(t.vals,ci.cut,type="l",lty=2)

#  Calculate Wald confidence interval.
that=tt0
theta1=theta10
theta2=theta20
I.mat=matrix(0,3,3)
I.mat[1,1]=1/(theta1*(1-theta1))
I.mat[2,2]=1/(theta2*(1-theta2))
I.mat[3,3]=(1-(1-theta1)*(1-theta2))/((1-theta1)*(1-theta2))
I.mat[1,3]=1/(1-theta1)
I.mat[3,1]=I.mat[1,3]
I.mat[2,3]=1/(1-theta2)
I.mat[3,2]=I.mat[2,3]
Vhat=solve(I.mat)
Vhat
me=1.96*sqrt(that*Vhat[3,3])
c((that-me),that,(that+me))