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


## -----------------------------------------------------------------------------
library(modelskill)

set.seed(123)

n <- 500

# Predictive means
pred <- seq(0, 10, length.out = n)

# True predictive standard deviation
predictive_sd <- rep(1, n)

# Independent observations generated from the predictive distributions
obs <- stats::rnorm(
  n,
  mean = pred,
  sd = predictive_sd
)


## -----------------------------------------------------------------------------
lower95 <- pred + stats::qnorm(0.025) * predictive_sd
upper95 <- pred + stats::qnorm(0.975) * predictive_sd


## -----------------------------------------------------------------------------
picp(obs, lower95, upper95)


## -----------------------------------------------------------------------------
coverage_error(
  obs,
  lower95,
  upper95,
  level = 0.95
)


## -----------------------------------------------------------------------------
interval_width(
  obs,
  lower95,
  upper95
)


## -----------------------------------------------------------------------------
interval_score(
  obs,
  lower95,
  upper95,
  level = 0.95
)


## -----------------------------------------------------------------------------
uncertainty_metrics(
  obs,
  lower = lower95,
  upper = upper95,
  level = 0.95
)

# Equivalent direct interface under the normal assumption:
uncertainty_metrics(obs, pred = pred, predictive_sd = predictive_sd, level = 0.95)


## -----------------------------------------------------------------------------
sd_too_small <- rep(0.5, n)

lower95_narrow <- pred + stats::qnorm(0.025) * sd_too_small
upper95_narrow <- pred + stats::qnorm(0.975) * sd_too_small

uncertainty_metrics(
  obs,
  lower = lower95_narrow,
  upper = upper95_narrow,
  level = 0.95
)


## -----------------------------------------------------------------------------
uncertainty_metrics(
  obs,
  lower = lower95,
  upper = upper95,
  level = 0.95
)


## -----------------------------------------------------------------------------
gg_coverage(
  obs,
  pred = pred,
  predictive_sd = predictive_sd
)


## -----------------------------------------------------------------------------
gg_coverage(
  obs,
  pred = pred,
  predictive_sd = sd_too_small
)


## -----------------------------------------------------------------------------
accuracy_plot_metrics(
  obs,
  pred = pred,
  predictive_sd = predictive_sd
)

## -----------------------------------------------------------------------------
accuracy_plot_metrics(
  obs,
  pred = pred,
  predictive_sd = sd_too_small
)

## -----------------------------------------------------------------------------
biased_pred <- pred + 0.5
biased_sd <- rep(1.118, n)

lower90_biased <-
  biased_pred + stats::qnorm(0.05) * biased_sd

upper90_biased <-
  biased_pred + stats::qnorm(0.95) * biased_sd

uncertainty_metrics(
  obs,
  lower = lower90_biased,
  upper = upper90_biased,
  level = 0.90
)


## -----------------------------------------------------------------------------
q_levels <- seq(0.05, 0.95, by = 0.05)

qhat <- vapply(
  q_levels,
  function(p) {
    pred + stats::qnorm(p) * predictive_sd
  },
  numeric(n)
)

qcp(
  obs,
  quantiles = qhat,
  levels = q_levels
)

## -----------------------------------------------------------------------------
gg_qcp(
  obs,
  quantiles = qhat,
  levels = q_levels
)

## -----------------------------------------------------------------------------
gg_qcp(
  obs,
  pred = pred,
  predictive_sd = predictive_sd
)

## -----------------------------------------------------------------------------
qhat_biased <- vapply(
  q_levels,
  function(p) {
    biased_pred + stats::qnorm(p) * biased_sd
  },
  numeric(n)
)

qcp(
  obs,
  quantiles = qhat_biased,
  levels = q_levels
)

gg_qcp(
  obs,
  quantiles = qhat_biased,
  levels = q_levels
)

## -----------------------------------------------------------------------------
gg_qcp(
  obs,
  pred = biased_pred,
  predictive_sd = biased_sd
)

## -----------------------------------------------------------------------------
pit_normal <- pit(
  obs = obs,
  pred = pred,
  predictive_sd = predictive_sd
)

gg_pit(pit_normal)

## -----------------------------------------------------------------------------
pit_narrow <- pit(
  obs = obs,
  pred = pred,
  predictive_sd = sd_too_small
)

gg_pit(pit_narrow)

## -----------------------------------------------------------------------------
pit_biased <- pit(
  obs = obs,
  pred = biased_pred,
  predictive_sd = biased_sd
)

gg_pit(pit_biased)

## -----------------------------------------------------------------------------
cdf_at_obs <- stats::pnorm(
  obs,
  mean = pred,
  sd = predictive_sd
)

pit_from_cdf <- pit(
  cdf_at_obs = cdf_at_obs
)

gg_pit(pit_from_cdf)

## -----------------------------------------------------------------------------
crps(
  obs,
  pred = pred,
  predictive_sd = predictive_sd
)

## -----------------------------------------------------------------------------
crps(
  obs,
  pred = pred,
  predictive_sd = predictive_sd
)

crps(
  obs,
  pred = pred,
  predictive_sd = sd_too_small
)

## -----------------------------------------------------------------------------
median_crps(
  obs,
  pred = pred,
  predictive_sd = predictive_sd
)

## -----------------------------------------------------------------------------
density_at_obs <- stats::dnorm(
  obs,
  mean = pred,
  sd = predictive_sd
)

log_score(
  obs,
  density_at_obs
)

## -----------------------------------------------------------------------------
set.seed(456)

n_draws <- 200

predictive_samples <- sapply(
  seq_len(n_draws),
  function(j) {
    stats::rnorm(
      n,
      mean = pred,
      sd = predictive_sd
    )
  }
)

dim(predictive_samples)


## -----------------------------------------------------------------------------
uncertainty_metrics(
  obs,
  distribution = predictive_samples,
  level = 0.95
)

# Individual interval statistics accept the same input.
picp(obs, distribution = predictive_samples, level = 0.95)
interval_width(obs, distribution = predictive_samples, level = 0.95)
interval_score(obs, distribution = predictive_samples, level = 0.95)


## -----------------------------------------------------------------------------
interval_levels <- seq(0.10, 0.90, by = 0.10)

gg_coverage(
  obs,
  distribution = predictive_samples,
  levels = interval_levels
)


## -----------------------------------------------------------------------------
accuracy_plot_metrics(
  obs,
  distribution = predictive_samples,
  levels = interval_levels
)


## -----------------------------------------------------------------------------
qcp(
  obs,
  distribution = predictive_samples,
  levels = q_levels
)


## -----------------------------------------------------------------------------
gg_qcp(
  obs,
  distribution = predictive_samples,
  levels = q_levels
)


## -----------------------------------------------------------------------------
pit_samples <- pit(obs = obs, distribution = predictive_samples)

gg_pit(pit_samples)


## -----------------------------------------------------------------------------
crps(
  obs,
  distribution = predictive_samples
)


## -----------------------------------------------------------------------------
median_crps(
  obs,
  distribution = predictive_samples
)


## -----------------------------------------------------------------------------
crps_decomposition(
  obs,
  distribution = predictive_samples
)


## -----------------------------------------------------------------------------
density_at_obs <- stats::dnorm(
  obs,
  mean = pred,
  sd = predictive_sd
)

log_score(obs, density_at_obs)