## ----setup, include=FALSE-------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
suppressPackageStartupMessages(library(fitPS))

## ----roux-data------------------------------------------------------
data("Psurveys")
roux = Psurveys$roux
roux

## ----mle-fits-------------------------------------------------------
mleFits = list(
  Zeta = fit(roux, model = zetaModel()),
  ZIZ = fit(roux, model = zizModel()),
  Logarithmic = fit(roux, model = logarithmicModel())
)

## ----mle-comparison-------------------------------------------------
mleComparison = do.call(
  rbind,
  lapply(names(mleFits), function(modelName) {
    fittedModel = mleFits[[modelName]]
    logLikelihood = logLik(fittedModel)

    data.frame(
      model = modelName,
      parameters = attr(logLikelihood, "df"),
      method = "MLE",
      engine = NA_character_,
      logLik = as.numeric(logLikelihood),
      AIC = AIC(fittedModel),
      BIC = BIC(fittedModel),
      DIC = NA_real_,
      check.names = FALSE
    )
  })
)

knitr::kable(mleComparison, digits = 2)

aicBest = mleComparison$model[which.min(mleComparison$AIC)]
bicBest = mleComparison$model[which.min(mleComparison$BIC)]

## ----bayesian-fits--------------------------------------------------
bayesFits = list(
  Zeta = fit(
    roux,
    model = zetaModel(),
    method = "bayes",
    bayesOptions = list(posteriorMethod = "numerical")
  ),
  ZIZ = fit(
    roux,
    model = zizModel(),
    method = "bayes",
    bayesOptions = list(posteriorMethod = "numerical")
  ),
  Logarithmic = fit(
    roux,
    model = logarithmicModel(),
    method = "bayes",
    bayesOptions = list(posteriorMethod = "numerical")
  )
)

## ----combined-comparison--------------------------------------------
bayesComparison = do.call(
  rbind,
  lapply(names(bayesFits), function(modelName) {
    data.frame(
      model = modelName,
      parameters = mleComparison$parameters[mleComparison$model == modelName],
      method = "Bayes",
      engine = "numerical",
      logLik = NA_real_,
      AIC = NA_real_,
      BIC = NA_real_,
      DIC = as.numeric(DIC(bayesFits[[modelName]])),
      check.names = FALSE
    )
  })
)

comparison = rbind(mleComparison, bayesComparison)
dicBest = bayesComparison$model[which.min(bayesComparison$DIC)]
knitr::kable(comparison, digits = 2)

## ----probability-table----------------------------------------------
maxTerm = max(roux$data$n) + 2L
terms = 0:maxTerm

observed = numeric(length(terms))
observed[match(roux$data$n, terms)] = roux$data$rn / sum(roux$data$rn)

probabilityComparison = data.frame(
  term = terms,
  Observed = observed,
  Zeta = as.numeric(predict(mleFits$Zeta, newdata = terms)),
  ZIZ = as.numeric(predict(mleFits$ZIZ, newdata = terms)),
  Logarithmic = as.numeric(predict(mleFits$Logarithmic, newdata = terms)),
  check.names = FALSE
)

knitr::kable(probabilityComparison, digits = 4)

## ----probability-plot, echo=FALSE-----------------------------------
plot(
  probabilityComparison$term,
  probabilityComparison$Observed,
  type = "h",
  lwd = 3,
  xlab = "Number of glass sources, k",
  ylab = "Probability",
  ylim = range(probabilityComparison[, -1])
)
points(
  probabilityComparison$term,
  probabilityComparison$Observed,
  pch = 16
)
matlines(
  probabilityComparison$term,
  as.matrix(probabilityComparison[, c("Zeta", "ZIZ", "Logarithmic")]),
  lty = 1:3,
  lwd = 2
)
legend(
  "topright",
  legend = c("Observed", "Zeta", "ZIZ", "Logarithmic"),
  pch = c(16, NA, NA, NA),
  lty = c(NA, 1:3),
  lwd = c(NA, 2, 2, 2),
  bty = "n"
)

## ----posterior-predictive-probabilities-----------------------------
predict(
  bayesFits$ZIZ,
  newdata = terms,
  type = "posteriorMean",
  interval = "credible"
)