library(RTMB)
dat <- read.table("bh.dat", header=TRUE)
par <- list(loga=0, logb=0, logsigma=0)

nll<-function(par){
  getAll(par, dat)
  sigma <- exp(logsigma)
  pred <- loga+log(ssb)-log(1+exp(logb)*ssb)
  nll <- -sum(dnorm(logR,pred,sigma,TRUE))
  return(nll)
}

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

ls.pro<-TMB:::tmbprofile(obj, "logsigma")
plot(ls.pro)
confint(ls.pro)

density.tmbprofile<-function(x,...){
  xx<-x[,1]
  yy<-exp(-x[,2])
  A<-sum(diff(xx) * (yy[-length(yy)]+yy[-1]))/2
  return(list(x=xx,y=yy/A, data.name=names(x)[1]))
}
pro <- TMB:::tmbprofile(obj,"loga",ytol=10)
sdrep<-summary(TMB:::sdreport(obj))
plot(pro[,1],dnorm(pro[,1],sdrep["loga",1], sdrep["loga",2]), 
     type="l", lwd=2, xlab="loga", ylab="")
lines(density(pro), lwd=2, col="red")




#