## ----quick-setup, include=FALSE-----------------------------
library(IJSE)
options(width = 62)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## ----quick-data---------------------------------------------
set.seed(431)
y <- rnorm(80, mean = 0.5, sd = 2)
n <- length(y)
v <- 1 / (n + 1 / 4)
m <- v * sum(y)
mu <- rnorm(20000, mean = m, sd = sqrt(v))
draws <- cbind(mu = mu)
rownames(draws) <- paste0("s", seq_along(mu))
ll <- vapply(y, function(yi) dnorm(yi, mu, log = TRUE),
             numeric(length(mu)))
rownames(ll) <- rownames(draws)
c(draw_rows = nrow(draws), likelihood_rows = nrow(ll),
  observations = ncol(ll))

## ----quick-fit----------------------------------------------
result <- ijse_draws(draws, ll, diagnostics = "none")
data.frame(Estimate = unname(coef(result)),
  Posterior_SD = unname(result$posterior_sd),
  IJ_SE = unname(result$se))

## ----quick-intervals----------------------------------------
confint(result, level = 0.95)
vcov(result)
vcov(result, type = "posterior")

## ----quick-oracle-------------------------------------------
exact_se <- v * sqrt(n * var(y))
data.frame(Analytic_IJ = exact_se,
           Draw_based_IJ = unname(result$se))

## ----quick-identity-----------------------------------------
stopifnot(isTRUE(all.equal(
  unname(result$se) / exact_se, var(mu) / v,
  tolerance = 1e-10)))

## ----quick-transform----------------------------------------
nonlinear <- ijse_draws(draws, ll,
  transform = function(d) {
    cbind(mu = d[, "mu"], mu_squared = d[, "mu"]^2)
  }, diagnostics = "none")
coef(nonlinear)
c(Exact_E_mu_squared = m^2 + v, Squared_E_mu = m^2)
vcov(nonlinear)

## ----quick-contrast-----------------------------------------
contrast <- ijse_draws(draws, ll,
  transform = function(d) d[, "mu"]^2 - 2 * d[, "mu"],
  diagnostics = "none")
a <- c(mu = -2, mu_squared = 1)
contrast_se <- sqrt(drop(t(a) %*% vcov(nonlinear) %*% a))
stopifnot(isTRUE(all.equal(
  unname(contrast$se), contrast_se)))

## ----quick-stream-------------------------------------------
read_likelihood <- function(ii) {
  block <- vapply(y[ii],
    function(yi) dnorm(yi, mu, log = TRUE),
    numeric(length(mu)))
  rownames(block) <- rownames(draws)
  block
}
streamed <- ijse_draws(draws, read_likelihood,
  n_obs = n, block_size = 7, diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(streamed), vcov(result))))

## ----quick-record-------------------------------------------
record <- list(
  estimate = coef(result),
  ij_covariance = vcov(result),
  posterior_covariance = vcov(result, type = "posterior"),
  interval = confint(result),
  draws = result$n_draws, observations = result$n_obs,
  units = result$n_units, clustered = result$clustered,
  version = as.character(packageVersion("IJSE")))
stopifnot(record$units == n, !record$clustered)