dat <- list()
#dat$X <- rmultinom(100, 15, c(.1,.5,.2,.05,.15))
dat$X <- cbind(rmultinom(50, 15, c(.1,.5,.2,.05,.15)), rmultinom(50, 15, rev(c(.1,.5,.2,.05,.15))))
library(RTMB)

par<-list(alpha=numeric(nrow(dat$X)-1))

trans <- function (alpha){
 expa <- exp(alpha)
 c(expa , 1+numeric(1))/(1+sum(expa))
}

nll <- function (par){
  getAll(par,dat)
  X<-OBS(X)  
  p <- trans(alpha)
  ret <- 0
  for(i in 1:ncol(X)){
    ret <- ret - dmultinom(X[,i], prob=p, log=TRUE)
  }
  ret
}

obj <- MakeADFun(nll,par)
opt <- nlminb (obj$par, obj$fn, obj$gr)

r<-oneStepPredict(obj, method="oneStepGeneric", discrete=TRUE, discreteSupport=0:15)$residual

stockassessment:::plotby(as.vector(col(dat$X)), as.vector(row(dat$X)), r); abline(h=nrow(dat$X), lwd=50, col="hotpink")