## ----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

## ----fit-mle--------------------------------------------------------
mleFit = fit(
  roux,
  model = zizModel()
)

## ----fit-bootstrap, eval=FALSE--------------------------------------
# bootFit = fit(
#   roux,
#   model = zizModel(),
#   method = "bootstrap",
#   B = 1000,
#   seed = 1234,
#   silent = TRUE
# )

## ----fit-bayes------------------------------------------------------
bayesFit = fit(
  roux,
  model = zizModel(),
  method = "bayes"
)

## ----fit-bayesian-bootstrap, eval=FALSE-----------------------------
# bayesBoot = fit(
#   roux,
#   model = zizModel(),
#   method = "bayesianBootstrap",
#   B = 1000,
#   seed = 1234
# )

## ----load-bootstrap-results, include=FALSE--------------------------
uncertaintyCache = readRDS(file.path("figures", "uncertainty-bootstrap-cache.rds"))
bootFit = uncertaintyCache$bootFit
bayesBoot = uncertaintyCache$bayesBoot

## ----mle-regions, results='hide'------------------------------------
mleRegions = confint(
  mleFit,
  level = c(0.80, 0.95)
)

## ----mle-region-plot, fig.cap="Profile-likelihood confidence regions for the ZIZ inflation and shape parameters."----
plotUncertainty(mleFit, level = c(0.80, 0.95))

## ----bayes-regions, results='hide'----------------------------------
bayesRegions = credint(
  bayesFit,
  level = c(0.80, 0.95),
  silent = TRUE
)

## ----uncertainty-region-comparison-code, eval=FALSE-----------------
# plotUncertainty(mleFit, level = c(0.80, 0.95))
# plotUncertainty(bootFit, level = c(0.80, 0.95))
# plotUncertainty(bayesFit, level = c(0.80, 0.95))
# plotUncertainty(bayesBoot, level = c(0.80, 0.95))

## ----uncertainty-region-comparison, echo=FALSE, out.width="100%", fig.cap="ZIZ parameter uncertainty under four inferential approaches. Top left: profile-likelihood confidence regions. Top right: ordinary-bootstrap regions with bootstrap realizations. Bottom left: parametric-Bayesian credible regions. Bottom right: Rubin Bayesian-Bootstrap regions with weighted-fit realizations."----
knitr::include_graphics(file.path("figures", "uncertainty-comparison.png"))

## ----mle-probabilities----------------------------------------------
fitted(mleFit, n = 6, type = "plugIn")

## ----bootstrap-probabilities----------------------------------------
bootstrapProbs(bootFit, n = 6)

## ----posterior-probabilities----------------------------------------
posteriorProbs(bayesFit, n = 6)

## ----bayesian-bootstrap-probabilities-------------------------------
bayesBoot$probabilities[1:6, ]

## ----compare-four-probability-estimates-----------------------------
probabilityComparison = data.frame(
  term = bootstrapProbs(bootFit, n = 6)$term,
  mlePlugIn = unname(fitted(mleFit, n = 6, type = "plugIn")),
  bootstrapMean = bootstrapProbs(bootFit, n = 6)$estimate,
  posteriorMean = posteriorProbs(bayesFit, n = 6)$estimate,
  bayesianBootstrapMean = bayesBoot$probabilities$estimate[1:6]
)

knitr::kable(
  probabilityComparison,
  digits = 5,
  col.names = c(
    "Term",
    "MLE plug-in",
    "Bootstrap mean",
    "Posterior mean",
    "Bayesian Bootstrap mean"
  )
)

## ----probability-uncertainty-comparison, echo=FALSE, fig.height=5, fig.cap="Uncertainty in the first six fitted P probabilities. The MLE plug-in value is a point estimate; the other three methods show their stored 95% uncertainty intervals."----
bootProb = bootstrapProbs(bootFit, n = 6)
posteriorProb = posteriorProbs(bayesFit, n = 6)
bayesBootProb = bayesBoot$probabilities[1:6, ]

terms = bootProb$term
methodNames = c(
  "MLE plug-in",
  "Ordinary bootstrap",
  "Parametric Bayes",
  "Bayesian Bootstrap"
)
methodOffset = c(-0.27, -0.09, 0.09, 0.27)
termPosition = seq_along(terms)

allLower = c(
  fitted(mleFit, n = 6, type = "plugIn"),
  bootProb$lower,
  posteriorProb$lower,
  bayesBootProb$lower
)
allUpper = c(
  fitted(mleFit, n = 6, type = "plugIn"),
  bootProb$upper,
  posteriorProb$upper,
  bayesBootProb$upper
)

plot(
  NA,
  xlim = range(c(allLower, allUpper), finite = TRUE),
  ylim = c(0.5, length(terms) + 0.5),
  yaxt = "n",
  ylab = "",
  xlab = "Fitted probability"
)
axis(2, at = termPosition, labels = terms, las = 1)

mleValues = unname(fitted(mleFit, n = 6, type = "plugIn"))
points(mleValues, termPosition + methodOffset[1], pch = 16)

intervalSets = list(bootProb, posteriorProb, bayesBootProb)
for (methodIndex in seq_along(intervalSets)) {
  intervalData = intervalSets[[methodIndex]]
  y = termPosition + methodOffset[methodIndex + 1L]
  segments(intervalData$lower, y, intervalData$upper, y)
  points(intervalData$estimate, y, pch = methodIndex)
}

legend(
  "topright",
  legend = methodNames,
  pch = c(16, 1, 2, 3),
  bty = "n",
  cex = 0.8
)

## ----probability-prediction-interface-------------------------------
predict(mleFit, newdata = 0:5, type = "plugIn")
predict(
  bootFit,
  newdata = 0:5,
  type = "bootstrapMean",
  interval = "percentile"
)
predict(
  bayesFit,
  newdata = 0:5,
  type = "posteriorMean",
  interval = "credible"
)

## ----practical-inflation--------------------------------------------
posteriorInflation(bayesFit, epsilon = 0.01)