mos.R
In mos: Simulation and Moment Computation for Order Statistics

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = TRUE,
  results = "markup",
  message = FALSE,
  warning = FALSE
)
library(mos)

## -----------------------------------------------------------------------------
ros(5, r = 2, n = 10, dist = "norm", mean = 0, sd = 1)

## -----------------------------------------------------------------------------
qgumbel <- function(p, mu = 0, beta = 1) mu - beta * log(-log(p))
ros(5, r = 3, n = 15, qf = "qgumbel", mu = 0, beta = 1)

## -----------------------------------------------------------------------------
ros(3, 1:5, 5, dist = "unif", min = 1, max = 5)

## -----------------------------------------------------------------------------
mo_unif(r = 2, n = 5, k = 1)  # Mean of 2nd order statistic

## -----------------------------------------------------------------------------
mo_exp(r = 3, n = 6, k = 2, mu = 0, sigma = 1)

## -----------------------------------------------------------------------------
mo_norm(r = 2, n = 10, k = 1)

## -----------------------------------------------------------------------------
mo_beta(r = 2, n = 5, k = 2, a = 2, b = 3)

## -----------------------------------------------------------------------------
rcens(n = 10, dist = "exp", type = "I", cens.time = 2, rate = 1)

## -----------------------------------------------------------------------------
rcens(n = 10, r = 5, dist = "norm", type = "II", mean = 0, sd = 1)

## -----------------------------------------------------------------------------
rpcens2(n = 10, R = c(2, 1, 2, 0, 0), dist = "exp", rate = 1)

## -----------------------------------------------------------------------------
rkrec(size = 5, k = 2, record = "upper", dist = "norm", mean = 0, sd = 1)

## -----------------------------------------------------------------------------
rkrec(size = 5, k = 3, record = "lower", dist = "exp", rate = 1)

## ----fig.width = 7.5, fig.height = 6.5----------------------------------------
if (requireNamespace("moments", quietly = TRUE)) {
  library(moments)
  set.seed(123)
  x <- ros(1e4, r = 1:25, n = 25, dist = "unif")
  observed <- colMeans(x^2)
  expected <- mo_unif(r = 1:25, n = 25, k = 2, a = 0, b = 1)
  oldpar <- par(mfrow = c(2, 2))
  on.exit(par(oldpar))
  plot(observed, expected, main = "2nd Moment of U(0,1)", xlab = "Observed", ylab = "Expected")
  abline(0, 1, col = 2)

  var_obs <- apply(x, 2, var)
  var_expc <- varOS(1:25, 25, dist = "unif", a = 0, b = 1)
  plot(var_obs, var_expc, main = "Variance of Order Statistics from U(0,1)",
       xlab = "Observed", ylab = "Expected")
  abline(0, 1, col = 2)

  skew_obs <- apply(x, 2, moments::skewness)
  skew_expc <- skewOS(1:25, 25, dist = "unif", a = 0, b = 1)
  plot(skew_obs, skew_expc, main = "Skewness of Order Statistics from U(0,1)",
       xlab = "Observed", ylab = "Expected")
  abline(0, 1, col = 2)

  kurt_obs <- apply(x, 2, moments::kurtosis)
  kurt_expc <- kurtOS(1:25, 25, dist = "unif", a = 0, b = 1)
  plot(kurt_obs, kurt_expc, main = "Kurtosis of Order Statistics from U(0,1)",
       xlab = "Observed", ylab = "Expected")
  abline(0, 1, col = 2)
}