library(bbmle)
library(numDeriv)
library(MCMCpack)
library(coda)

sessionInfo()

ReedfrogSizepred <- read.table("../DATA/tadpole.dat",header=TRUE)

source("tadpole_R_funs.R")

t.fit <- system.time(tfit <- tadpole_R_fit())
t.profile <- system.time(tprof <- tadpole_R_profile(tfit))
t.mcmc <- system.time(tmcmc <- tadpole_R_mcmc(tfit))

times <- c(fit=t.fit["elapsed"],profile=t.profile["elapsed"],mcmc=t.mcmc["elapsed"])
fits <- list(fit=tfit,profile=tprof,mcmc=tmcmc,times=times)
save("fits",file="fit_fits.RData")

## FIXME: can replace with @details$maxgrad, @details$eratio
##  when bbmle 1.0.1 is available from R-forge/CRAN ...

## compute termination info
maxgrad <-
  max(abs(grad(function(x)
       do.call(tfit@minuslogl,as.list(x)),
       coef(tfit))))

ev <- eigen(tfit@details$hessian)$value
eratio <- min(ev)/max(ev)

results <- list(obj=-logLik(tfit),  ## objective function/Neg log likelihood (NULL for MCMC)
                coef=coef(tfit),  ## coefficients
                sd=sqrt(diag(vcov(tfit))),
                vcov=vcov(tfit),
                confint.quad=confint(tfit,method="quad"),
                confint.profile=confint(tprof),
                confint.mcmc=t(apply(tmcmc,2,quantile,c(0.025,0.975))),
                time=times,
                terminfo=c(maxgrad=maxgrad,eratio=eratio), ## MLE termination info
                convinfo=list(effsize=effectiveSize(tmcmc),
                  geweke=geweke.diag(tmcmc))## MCMC convergence info
                )

save("results",file="fitx.RData")