ovenCounts.MLE = function() {
    
# define log likelihood function
logLike = function(yMat, nrepls, lambda, p, NinftyMaximum=5000) {
  nsites = length(nrepls)
  logLike = rep(NA, nsites)
  for (i in 1:nsites) {
    y = yMat[i, 1:nrepls[i]]
    Nmin = max(y)
    
    # compute normalizing constant for Poisson
    Ninfty = Nmin
    deltaNinfty = 10
    denom = 0
    while (abs(1-denom)>=0.001 & Ninfty<NinftyMaximum) {
      Ninfty = Ninfty + deltaNinfty
      supportOfN = 0:Ninfty
      denom = sum(exp(dpois(supportOfN, lambda=lambda[i], log=T)))
    }
    if (Ninfty >= NinftyMaximum) {
      cat(Ninfty, " (= upper bound for support of Poisson) exceeds maximum allowable value \n")
    }
    
    # compute marginal likelihood for site
    siteSum = 0
    if (Nmin==0) {
      siteSum = siteSum + exp(-lambda[i])
      Nmin=1
    }
    N = Nmin:Ninfty
    logSum =  dpois(N, lambda=lambda[i], log=T)
    for (j in 1:nrepls[i]) {
      logSum = logSum + dbinom(y[j], size=N, prob=p[i], log=T)
    }
    siteSum = siteSum + sum(exp(logSum))
    
    logLike[i] = log(siteSum)
  }
  
  sum(logLike)
}


# define negative log likeihood function
  negLL = function(param) {
    pVec = rep(expit(param[1]), nsites)
    lambdaVec = rep(exp(param[2]), nsites)
    (-1)*logLike(yMat, nrepls, lambdaVec, pVec)
  }
  

# read data from file and pre-process for analysis
  yMat = as.matrix(read.table('oven.txt'))

  nsites = dim(yMat)[1]
  nrepls= rep(dim(yMat)[2], nsites)
  yMax = rep(NA, nsites)
  for (i in 1:nsites) {
    yMax[i] = max(yMat[i, 1:nrepls[i]])
  }
  pGuess = 0.5
  lambdaGuess = mean(yMax)/pGuess


# fit model to data
  fit = optim(par=c(logit(pGuess), log(lambdaGuess)), fn=negLL, method='BFGS', hessian=T)
  
  list(logLikelihood=-fit$value, mle=fit$par, covMat=chol2inv(chol(fit$hessian)))
}