---
title: "Working with Order Statistics in the mos Package"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Working with Order Statistics in the mos Package}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
echo = TRUE,
results = "markup",
message = FALSE,
warning = FALSE
)
library(mos)
```
## Introduction
The `mos` package provides a suite of tools to work with order statistics, including:
* Generating random values of order statistics and censored/record samples,
* Analytical computation of moments (mean, variance, skewness, kurtosis),
* Simulation-based estimation of moments for some of continuous distributions.
These tools are useful for researchers in statistics, reliability, and simulation studies.
---
## Generating Order Statistics
You can simulate order statistics from continuous distributions using the `ros()` function. This function accepts either:
* A named distribution available in base R (e.g., "norm", "exp"), or
* A user-defined quantile function using the `qf` argument.
### Example 1: Using a built-in distribution
```{r}
ros(5, r = 2, n = 10, dist = "norm", mean = 0, sd = 1)
```
### Example 2: Using a custom quantile function (e.g., Gumbel)
```{r}
qgumbel <- function(p, mu = 0, beta = 1) mu - beta * log(-log(p))
ros(5, r = 3, n = 15, qf = "qgumbel", mu = 0, beta = 1)
```
### Example 3: Generate multiple order statistics from the uniform distribution
```{r}
ros(3, 1:5, 5, dist = "unif", min = 1, max = 5)
```
---
## Moments of Order Statistics
The package provides two categories of moment computation:
* **Analytical formulas** (closed-form), where available
* **Simulation-based estimation** otherwise
### Distributions with analytical support
The following distributions have closed-form expressions for moments:
* Uniform
* Exponential
* Weibull
* Symmetric Triangular
* Topp-Leon
* Complementary Beta (for integer `b` and `k <= 2`)
### Distributions handled via simulation
* Normal
* Student's t
* Gamma
* Beta
* Pareto
* Kumaraswamy
* Complementary Beta (for non-integer `b` or `k > 2`)
In all simulation-based cases, the `ros()` function is used internally to generate the order statistics.
### Example: Uniform Distribution (Analytical)
```{r}
mo_unif(r = 2, n = 5, k = 1) # Mean of 2nd order statistic
```
### Example: Exponential Distribution (Analytical)
```{r}
mo_exp(r = 3, n = 6, k = 2, mu = 0, sigma = 1)
```
### Example: Normal Distribution (Simulated)
```{r}
mo_norm(r = 2, n = 10, k = 1)
```
### Example: Beta Distribution (Simulated)
```{r}
mo_beta(r = 2, n = 5, k = 2, a = 2, b = 3)
```
In addition to `mo_*()` functions, the package includes:
* `varOS()` to compute variance
* `skewOS()` to compute skewness
* `kurtOS()` to compute kurtosis
---
## Censored Samples
The `mos` package supports simulation of censored data using several standard censoring schemes:
### Type I Censoring (Time-based)
Observations are censored beyond a fixed time limit.
```{r}
rcens(n = 10, dist = "exp", type = "I", cens.time = 2, rate = 1)
```
### Type II Censoring (Failure-based)
The smallest `r` values are observed (i.e., censoring after a fixed number of failures).
```{r}
rcens(n = 10, r = 5, dist = "norm", type = "II", mean = 0, sd = 1)
```
### Progressive Type-II Censoring
Implements progressive removal of remaining items based on a censoring scheme.
```{r}
rpcens2(n = 10, R = c(2, 1, 2, 0, 0), dist = "exp", rate = 1)
```
This is useful in reliability experiments where censoring occurs at each failure stage.
---
## Record Values
The `rkrec()` function generates upper and lower k-records from continuous distributions.
### Upper k-Records
```{r}
rkrec(size = 5, k = 2, record = "upper", dist = "norm", mean = 0, sd = 1)
```
### Lower k-Records
```{r}
rkrec(size = 5, k = 3, record = "lower", dist = "exp", rate = 1)
```
---
## An Application of the Package in Graphical Goodness-of-Fit
This section compares analytical and simulated values of the second moment, variance, skewness, and kurtosis of order statistics from the Uniform(0,1) distribution:
```{r, 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)
}
```
The plots show a strong alignment between the observed and expected values along the identity line. This suggests that the order statistics simulated from the Uniform(0,1) distribution closely follow their theoretical moments. Such alignment in graphical comparisons is often taken as visual evidence of a good fit between observed data and the theoretical model.
---
## Conclusion
The `mos` package integrates analytical and simulation-based methods for working with order statistics, censored data, and record values. It serves as a powerful and flexible toolkit for statisticians, particularly those involved in reliability analysis and the study of ordered data.