| Title: | Functions and Data Sets for "That's Weird: Anomaly Detection Using R" by Rob J Hyndman |
|---|---|
| Description: | All functions and data sets required for the examples in the book Hyndman (2026) "That's Weird: Anomaly Detection Using R" <https://OTexts.com/weird/>. All packages needed to run the examples are also loaded. |
| Authors: | Rob Hyndman [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-2140-5352>), Torben Tvedebrink [ctb], Posit Software, PBC [cph] |
| Maintainer: | Rob Hyndman <[email protected]> |
| License: | GPL-3 |
| Version: | 3.1.0.9000 |
| Built: | 2026-10-14 07:58:34 UTC |
| Source: | https://github.com/robjhyndman/weird |
Augment the data with information from an rrcov::Pca* object
(such as the output of rrcov::PcaHubert() or rrcov::PcaClassic()). The
returned tibble contains the principal component scores (.fittedPC1,
.fittedPC2, ...), the score distance (.sd) and the orthogonal distance
(.od) of each observation. The score distance measures how far an
observation lies from the centre within the projection subspace, while the
orthogonal distance measures how far it lies from the subspace. If data
is supplied, its columns are returned alongside these results.
## S3 method for class 'Pca' augment(x, data = NULL, ...)## S3 method for class 'Pca' augment(x, data = NULL, ...)
x |
An |
data |
The original data matrix or data frame used to compute the PCA. If supplied, its columns are bound to the left of the returned tibble. |
... |
Unused. |
A tibble::tibble() with one row per observation.
Rob J Hyndman
Y <- oldfaithful[, c("duration", "waiting")] pca <- rrcov::PcaHubert(as.matrix(Y), k = 1) broom::augment(pca, data = Y)Y <- oldfaithful[, c("duration", "waiting")] pca <- rrcov::PcaHubert(as.matrix(Y), k = 1) broom::augment(pca, data = Y)
Draw a two-dimensional projection of the scores with the
original variable axes overlaid as arrows (loadings), as in a biplot. Pass
object (the output of stats::prcomp() or an rrcov::Pca* function);
otherwise supply scores and loadings directly. All scores should be
centred about the origin. The arrows are stretched by a common factor so that
the longest arrow just reaches the edge of the point cloud, and only loadings
longer than label_threshold are labelled.
biplot_projection( object = NULL, scores = NULL, loadings = NULL, label_threshold = 0, arrow_colour = "#c14b14", ... )biplot_projection( object = NULL, scores = NULL, loadings = NULL, label_threshold = 0, arrow_colour = "#c14b14", ... )
object |
Optionally, the output of |
scores |
A matrix or data frame of scores centred about the origin,
with the first two columns used as the horizontal and vertical coordinates.
Ignored if |
loadings |
A matrix or data frame of loadings, with row names giving the
variable names and the first two columns used as the arrow directions.
Ignored if |
label_threshold |
Only loadings whose absolute length exceeds this
threshold are labelled. The default of |
arrow_colour |
Colour of the arrows and labels. |
... |
Additional arguments passed to |
A ggplot object.
Rob J Hyndman
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 9, https://OTexts.com/weird/highdim.html.
oldfaithful[, c("duration", "waiting")] |> prcomp(scale = TRUE) |> biplot_projection()oldfaithful[, c("duration", "waiting")] |> prcomp(scale = TRUE) |> biplot_projection()
A dataset containing career batting statistics for all international test players (men and women) up to 6 October 2025.
cricket_battingcricket_batting
A data frame with 3968 rows and 15 variables:
Player name in form of "initials surname"
Country played for
First year of test playing career
Last year of test playing career
Number of matches played
Number of innings batted
Number of times not out
Total runs scored
Highest score in an innings
Was highest score not out?
Batting average at end of career
Total number of 100s scored
Total number of 50s scored
Total number of 0s scored
"Men" or "Women"
Data frame
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 1.4, https://OTexts.com/weird/.
cricket_batting |> filter(Innings > 20) |> select(Player, Country, Matches, Runs, Average, Hundreds, Fifties, Ducks) |> arrange(desc(Average))cricket_batting |> filter(Innings > 20) |> select(Player, Country, Matches, Runs, Average, Hundreds, Fifties, Ducks) |> arrange(desc(Average))
Make a long-form data frame containing densities from a distributional object on a regular grid for plotting.
density_df(object, ngrid = NULL)density_df(object, ngrid = NULL)
object |
A distributional object |
ngrid |
Number of grid points in each dimension. Defaults to 501 for univariate distributions and 101 for bivariate distributions. |
Data frame with columns x, y (if bivariate), density, and distribution.
dist_kde(oldfaithful$duration) |> density_df()dist_kde(oldfaithful$duration) |> density_df()
Creates a distributional object using a kernel density estimate with a
Gaussian kernel obtained from the kde() function. The bandwidth
can be specified; otherwise the kde_bandwidth() function is used.
The cdf, quantiles and moments are consistent with the kde. Generating
random values from the kde is equivalent to a smoothed bootstrap.
dist_kde( y, h = NULL, H = NULL, method = c("robust", "normal", "plugin", "scv", "lookout"), ... )dist_kde( y, h = NULL, H = NULL, method = c("robust", "normal", "plugin", "scv", "lookout"), ... )
y |
Numerical vector or matrix of data, or a list of such objects. If a list is provided, then all objects should be of the same dimension. e.g., all vectors, or all matrices with the same number of columns. |
h |
Bandwidth for univariate distribution. Ignored if |
H |
Bandwidth matrix for multivariate distribution. If |
method |
The method of bandwidth estimation to use. See |
... |
Other arguments are passed to |
A distributional object of class dist_kde.
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 2.9 and 3.9, https://OTexts.com/weird/.
dist_kde(c(rnorm(200), rnorm(100, 5))) dist_kde(cbind(rnorm(200), rnorm(200, 5)), binned = TRUE)dist_kde(c(rnorm(200), rnorm(100, 5))) dist_kde(cbind(rnorm(200), rnorm(200, 5)), binned = TRUE)
Convert Gaussian mixture model to a distributional object
dist_mclust(object)dist_mclust(object)
object |
Object of class |
An object of class distributional
library(mclust, exclude = "count") # Univariate mixture density gmm <- Mclust(oldfaithful$duration) |> dist_mclust() gg_density(gmm) + geom_jitter(data = oldfaithful, aes(x=duration, y = -0.0002), width=0, height=0.0002, alpha = 0.1) + labs(x = "Eruption duration", y = "")library(mclust, exclude = "count") # Univariate mixture density gmm <- Mclust(oldfaithful$duration) |> dist_mclust() gg_density(gmm) + geom_jitter(data = oldfaithful, aes(x=duration, y = -0.0002), width=0, height=0.0002, alpha = 0.1) + labs(x = "Eruption duration", y = "")
All 1000 images of sneakers from the Fashion-MNIST test set, together with 10 images from other classes of varying similarity to sneakers: four ankle boots, three sandals, two bags and one pair of trousers. Each image is 28 x 28 pixels.
fashionfashion
A data frame with 1010 rows and 4 columns:
Index of the image in the Fashion-MNIST test set
Class of the image (a factor with the 10 Fashion-MNIST classes)
TRUE for the 10 images that are not sneakers
A 1010 x 784 integer matrix of pixel intensities, from 0
(white) to 255 (black), with one row per image. Each row contains the
28 x 28 pixels of an image, stored row by row, so
matrix(pixels[i, ], nrow = 28, byrow = TRUE) gives image i.
Data frame
Xiao, H., Rasul, K., & Vollgraf, R. (2017). Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms. https://github.com/zalandoresearch/fashion-mnist (MIT licence)
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 13, https://OTexts.com/weird/imagevideo.html.
fashion |> count(label)fashion |> count(label)
Hourly air quality measurements from 12 monitoring stations across Beijing, China, from 1 March 2013 to 28 February 2017. The data are downloaded and returned.
fetch_air_quality()fetch_air_quality()
A data frame with 420,768 rows and 17 columns:
Name of the monitoring station
Year of measurement
Month of measurement
Day of measurement
Hour of measurement (0–23)
Particulate matter with diameter less than 2.5 micrometers (micrograms per cubic meter)
Particulate matter with diameter less than 10 micrometers (micrograms per cubic meter)
Sulfur dioxide concentration (micrograms per cubic meter)
Nitrogen dioxide concentration (micrograms per cubic meter)
Carbon monoxide concentration (micrograms per cubic meter)
Ozone concentration (micrograms per cubic meter)
Temperature (degrees Celsius)
Atmospheric pressure (hPa)
Dew point temperature (degrees Celsius)
Rainfall (millimeters)
Wind direction
Wind speed (meters per second)
Data frame
Chen, S. (2017). Beijing Multi-Site Air Quality Dataset. UCI Machine Learning Repository. https://doi.org/10.24432/C5RK5G
Hyndman, R J (2026) That's weird: Anomaly detection using R, https://OTexts.com/weird/.
## Not run: air_quality <- fetch_air_quality() air_quality |> filter(station == "Aotizhongxin") |> ggplot(aes(x = temperature, y = pm2_5)) + geom_point(alpha = 0.1) ## End(Not run)## Not run: air_quality <- fetch_air_quality() air_quality |> filter(station == "Aotizhongxin") |> ggplot(aes(x = temperature, y = pm2_5)) + geom_point(alpha = 0.1) ## End(Not run)
The text of 29 Oz books: the 14 written by L. Frank Baum (1900–1920), and 15 of those written by Ruth Plumly Thompson (1921–1939). Each row is a paragraph. Front matter, chapter headings and titles, illustration captions and back matter have been removed, and curly quotes replaced by straight quotes. The data are downloaded and returned.
fetch_oz_books()fetch_oz_books()
A data frame with 31,258 rows and 7 columns:
Title of the book
Author of the book
Year of first publication
Project Gutenberg ebook number
Chapter number
Paragraph number within the book
Text of the paragraph
The Royal Book of Oz (1921) was credited to Baum, who had died in 1919, but was written by Thompson. The Magic of Oz (1919) and Glinda of Oz (1920) were published after Baum's death. Thompson's Oz books from 1931 to 1934 are not included as they are not available from Project Gutenberg.
Data frame
Project Gutenberg, https://www.gutenberg.org. Downloaded 9 October 2026.
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 10, https://OTexts.com/weird/.
## Not run: oz_books <- fetch_oz_books() oz_books |> group_by(title, author, year) |> summarise(paragraphs = n(), .groups = "drop") |> arrange(year) ## End(Not run)## Not run: oz_books <- fetch_oz_books() oz_books |> group_by(title, author, year) |> summarise(paragraphs = n(), .groups = "drop") |> arrange(year) ## End(Not run)
Greyscale frames from the PETS2006 video in the CDnet 2014 change detection benchmark. The video was taken by a fixed camera overlooking the concourse of a railway station. For most of the video, a man stands near the glass wall at the back of the concourse with a bag, while other people walk past. He then puts the bag down, and eventually walks away, leaving it behind. The original 720 x 576 colour frames have been converted to greyscale and reduced to 120 x 96 pixels by averaging over blocks of 6 x 6 pixels. The data are downloaded and returned.
fetch_pets2006()fetch_pets2006()
A data frame with 1200 rows and 2 columns:
Frame number
A 1200 x 11520 integer matrix of pixel intensities, from 0
(black) to 255 (white), with one row per frame. Each row contains the
96 x 120 pixels of a frame, stored column by column, so
matrix(pixels[i, ], nrow = 96) gives frame i as an image.
The benchmark uses frames 1–299 to initialise background models, and evaluates methods on frames 300–1200.
Data frame
Wang, Y., Jodoin, P.-M., Porikli, F., Konrad, J., Benezeth, Y., & Ishwar, P. (2014). CDnet 2014: An expanded change detection benchmark dataset. IEEE Conference on Computer Vision and Pattern Recognition Workshops, 393–400. http://changedetection.net
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 13, https://OTexts.com/weird/imagevideo.html.
## Not run: pets2006 <- fetch_pets2006() # Show frame 951 image(t(matrix(pets2006$pixels[951, ], nrow = 96))[, 96:1], col = grey.colors(256)) ## End(Not run)## Not run: pets2006 <- fetch_pets2006() # Show frame 951 image(t(matrix(pets2006$pixels[951, ], nrow = 96))[, 96:1], col = grey.colors(256)) ## End(Not run)
A data set containing data on wines from 42 countries, taken from Wine Enthusiast Magazine during the week of 15 June 2017. The data are downloaded and returned.
fetch_wine_reviews()fetch_wine_reviews()
A data frame with 111,592 rows and 9 columns:
Country of origin
State or province of origin
Region of origin
Name of vineyard that made the wine
Variety of grape
Year of wine
Points allocated by WineEnthusiast reviewer on a scale of 0-100
Price of a bottle of wine in $US
Text of the review
Data frame
https://www.kaggle.com/datasets/zynicide/wine-reviews
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 1.4, https://OTexts.com/weird/.
## Not run: wine_reviews <- fetch_wine_reviews() wine_reviews |> ggplot(aes(x = points, y = price)) + geom_jitter(height = 0, width = 0.2, alpha = 0.1) + scale_y_log10() ## End(Not run)## Not run: wine_reviews <- fetch_wine_reviews() wine_reviews |> ggplot(aes(x = points, y = price)) + geom_jitter(height = 0, width = 0.2, alpha = 0.1) + scale_y_log10() ## End(Not run)
A data set containing French mortality rates between the years 1816 and 1999, by age and sex.
fr_mortalityfr_mortality
A data frame with 31648 rows and 4 columns.
Data frame
Human Mortality Database https://www.mortality.org
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 1.4, https://OTexts.com/weird/.
fr_mortalityfr_mortality
Produces a bivariate bagplot. A bagplot is analogous to a univariate boxplot, except it is in two dimensions. Like a boxplot, it shows the median, a region containing 50% of the observations, a region showing the remaining observations other than outliers, and any outliers.
gg_bagplot(data, var1, var2, color = "#00659e", show_points = FALSE, ...)gg_bagplot(data, var1, var2, color = "#00659e", show_points = FALSE, ...)
data |
A data frame or matrix containing the data. |
var1 |
The name of the first variable to plot (a bare expression). |
var2 |
The name of the second variable to plot (a bare expression). |
color |
The base color to use for the median. Other colors are generated
as a mixture of |
show_points |
A logical argument indicating if a regular bagplot is required
( |
... |
Other arguments are passed to the |
A ggplot object showing a bagplot or scatterplot of the data.
Rob J Hyndman
Rousseeuw, P J, Ruts, I, & Tukey, J W (1999). The bagplot: A bivariate boxplot. The American Statistician, 52(4), 382–387.
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 5.6, https://OTexts.com/weird/.
gg_bagplot(n01, v1, v2) gg_bagplot(n01, v1, v2, show_points = TRUE)gg_bagplot(n01, v1, v2) gg_bagplot(n01, v1, v2, show_points = TRUE)
Produce ggplot of densities from distributional objects in 1 or 2 dimensions
gg_density( object, prob = seq(9)/10, hdr = NULL, show_points = FALSE, show_mode = FALSE, show_anomalies = FALSE, colors = c("#0072b2", "#D55E00", "#009E73", "#CC79A7", "#E69F00", "#56B4E9", "#F0E442", "#333333"), alpha = NULL, jitter = FALSE, ngrid = NULL )gg_density( object, prob = seq(9)/10, hdr = NULL, show_points = FALSE, show_mode = FALSE, show_anomalies = FALSE, colors = c("#0072b2", "#D55E00", "#009E73", "#CC79A7", "#E69F00", "#56B4E9", "#F0E442", "#333333"), alpha = NULL, jitter = FALSE, ngrid = NULL )
object |
distribution object from the distributional package or
|
prob |
Probability of the HDRs to be drawn. |
hdr |
Character string describing how the HDRs are to be shown. Options
are "none", "fill", "points" and "contours" (the latter only for bivariate plots).
If |
show_points |
If |
show_mode |
If |
show_anomalies |
If |
colors |
Color palette to use. If there are more than
|
alpha |
Transparency of points. Ignored if |
jitter |
For univariate distributions, when |
ngrid |
Number of grid points in each dimension, passed to
|
This function produces a ggplot of a density from a distributional object.
For univariate densities, it produces a line plot of the density function, with
an optional ribbon showing some highest density regions (HDRs) and/or the observations.
For bivariate densities, it produces an HDR contour plot of the density function, with
the observations optionally shown as points.
The mode can also be drawn as a point.
The combination of hdr = "fill", show_points = TRUE,
show_mode = TRUE, and prob = c(0.5, 0.99) is equivalent to showing
HDR boxplots.
A ggplot object.
Rob J Hyndman
# Univariate densities kde <- dist_kde(c(rnorm(500), rnorm(500, 4, 0.5))) gg_density(kde, hdr = "fill", prob = c(0.5, 0.95), color = "#c14b14", show_mode = TRUE, show_points = TRUE, jitter = TRUE ) c(dist_normal(), kde) |> gg_density(hdr = "fill", prob = c(0.5, 0.95)) # Bivariate density tibble(y1 = rnorm(5000), y2 = y1 + rnorm(5000)) |> dist_kde() |> gg_density(show_points = TRUE, alpha = 0.1, hdr = "fill")# Univariate densities kde <- dist_kde(c(rnorm(500), rnorm(500, 4, 0.5))) gg_density(kde, hdr = "fill", prob = c(0.5, 0.95), color = "#c14b14", show_mode = TRUE, show_points = TRUE, jitter = TRUE ) c(dist_normal(), kde) |> gg_density(hdr = "fill", prob = c(0.5, 0.95)) # Bivariate density tibble(y1 = rnorm(5000), y2 = y1 + rnorm(5000)) |> dist_kde() |> gg_density(show_points = TRUE, alpha = 0.1, hdr = "fill")
Produces a 1d or 2d box plot of HDR regions. The darker regions contain observations with higher probability, while the lighter regions contain points with lower probability. Observations outside the largest HDR are shown as individual points. Anomalies with leave-one-out surprisal probabilities less than 0.005 are optionally shown in black.
gg_hdrboxplot( data, var1, var2 = NULL, prob = c(0.5, 0.99), color = "#0072b2", show_points = FALSE, show_anomalies = TRUE, alpha = NULL, jitter = TRUE, ... )gg_hdrboxplot( data, var1, var2 = NULL, prob = c(0.5, 0.99), color = "#0072b2", show_points = FALSE, show_anomalies = TRUE, alpha = NULL, jitter = TRUE, ... )
data |
A data frame or matrix containing the data. |
var1 |
The name of the first variable to plot (a bare expression). |
var2 |
Optionally, the name of the second variable to plot (a bare expression). |
prob |
A numeric vector specifying the coverage probabilities for the HDRs. |
color |
The base color to use for the mode. Colors for the HDRs are generated by whitening this color. |
show_points |
A logical argument indicating if a regular HDR plot is required
( |
show_anomalies |
A logical argument indicating if the surprisal anomalies should be shown (in black). These are points with leave-one-out surprisal probability values less than 0.005 (using a GPD approximation), and which lie outside the 99% HDR region. |
alpha |
Transparency of points. Ignored if |
jitter |
A logical value indicating if the points should be vertically jittered for the 1d box plots to reduce overplotting. |
... |
Other arguments passed to |
The original HDR boxplot proposed by Hyndman (1996), can be produced
with show_anomalies = FALSE, jitter = FALSE, alpha = 1, and all other
arguments set to their defaults.
A ggplot object showing an HDR plot or scatterplot of the data.
Rob J Hyndman
Hyndman, R J (1996) Computing and Graphing Highest Density Regions, The American Statistician, 50(2), 120–126. https://robjhyndman.com/publications/hdr/
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 5.7, https://OTexts.com/weird/.
df <- data.frame(x = c(rnorm(1000), rnorm(1000, 5, 1), 10)) gg_hdrboxplot(df, x, show_anomalies = TRUE) cricket_batting |> filter(Innings > 20) |> gg_hdrboxplot(Average) oldfaithful |> gg_hdrboxplot(duration, waiting, show_points = TRUE)df <- data.frame(x = c(rnorm(1000), rnorm(1000, 5, 1), 10)) gg_hdrboxplot(df, x, show_anomalies = TRUE) cricket_batting |> filter(Innings > 20) |> gg_hdrboxplot(Average) oldfaithful |> gg_hdrboxplot(duration, waiting, show_points = TRUE)
Compute Global-Local Outlier Score from Hierarchies. This is based
on hierarchical clustering, using core distances to the k-th nearest neighbour. The resulting
outlier score is a measure of how anomalous each observation is.
The function uses dbscan::hdbscan to do the calculation.
glosh_scores(y, k = 10, ...)glosh_scores(y, k = 10, ...)
y |
Numerical matrix or vector of data |
k |
Number of neighbours to include, not counting the observation itself. Default: 10. |
... |
Additional arguments passed to |
Numerical vector containing GLOSH values
Rob J Hyndman
dbscan::glosh
y <- c(rnorm(49), 5) glosh_scores(y)y <- c(rnorm(49), 5) glosh_scores(y)
Grubbs' test (proposed in 1950) identifies possible anomalies in univariate data using z-scores assuming the data come from a normal distribution. Dixon's test (also from 1950) compares the difference in the largest two values to the range of the data. Critical values for Dixon's test have been computed using simulation with interpolation using a quadratic model on logit(alpha) and log(log(n)).
grubbs_anomalies(y, alpha = 0.05) dixon_anomalies(y, alpha = 0.05, two_sided = TRUE)grubbs_anomalies(y, alpha = 0.05) dixon_anomalies(y, alpha = 0.05, two_sided = TRUE)
y |
numerical vector of observations |
alpha |
size of the test. |
two_sided |
If |
Grubbs' test is based on z-scores, and a point is identified as an
anomaly when the associated absolute z-score is greater than a threshold value.
A vector of logical values is returned, where TRUE indicates an anomaly.
This version of Grubbs' test looks for outliers anywhere in the sample.
Grubbs' original test came in several variations which looked for one outlier,
or two outliers in one tail, or two outliers on opposite tails. These variations
are implemented in the grubbs.test function.
Dixon's test only considers the maximum (and possibly the minimum) as potential outliers.
A logical vector
Rob J Hyndman
Grubbs, F E (1950). Sample criteria for testing outlying observations. Annals of Mathematical Statistics, 21(1), 27–58.
Dixon, W J (1950). Analysis of extreme values. Annals of Mathematical Statistics, 21(4), 488–506.
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 4.4-4.5, https://OTexts.com/weird/.
x <- c(rnorm(1000), 5:10) tibble(x = x) |> filter(grubbs_anomalies(x)) tibble(x = x) |> filter(dixon_anomalies(x)) y <- c(rnorm(1000), 5) tibble(y = y) |> filter(grubbs_anomalies(y)) tibble(y = y) |> filter(dixon_anomalies(y))x <- c(rnorm(1000), 5:10) tibble(x = x) |> filter(grubbs_anomalies(x)) tibble(x = x) |> filter(dixon_anomalies(x)) y <- c(rnorm(1000), 5) tibble(y = y) |> filter(grubbs_anomalies(y)) tibble(y = y) |> filter(dixon_anomalies(y))
A data set containing gun ownership rates and homicide rates for 2017 for various countries around the world. The gun ownership rates are the number of guns owned by civilians per 100 people. The homicide rates are the number of homicides per 100,000 people where the weapon was a firearm.
gun_deathsgun_deaths
A data frame with 77 rows and 4 columns:
Country name
World region according to Our World in Data
Gun ownership rate (number of guns owned by civilians per 100 people)
Homicide rate (number of homicides per 100,000 people where the weapon was a firearm)
Data frame
World Population Review https://worldpopulationreview.com/country-rankings/gun-ownership-by-country and https://ourworldindata.org/grapher/homicide-rates-from-firearms
gun_deathsgun_deaths
The Hampel filter is designed to find anomalies in time series data using mean absolute deviations in the vicinity of each observation.
hampel_anomalies( y, bandwidth, k = 3, alpha = NULL, approximation = c("none", "gpd", "empirical") )hampel_anomalies( y, bandwidth, k = 3, alpha = NULL, approximation = c("none", "gpd", "empirical") )
y |
numeric vector containing time series |
bandwidth |
integer specifying half-width of the window around each observation. That is, each window contains observations from t-bandwidth, ..., t+bandwidth. Must be at least 1. |
k |
numeric number of standard deviations to declare an anomaly. Ignored if |
alpha |
numeric significance level for declaring an anomaly under a normal distribution.
If specified, |
approximation |
character string specifying the method to use for approximating the tail of the distribution of surprisal values. Options are "none" (no approximation), "gpd" (generalized Pareto distribution), or "empirical" (empirical distribution). |
First, a moving median is calculated using windows of size
2 * bandwidth + 1. Then the median absolute deviations from
this moving median are calculated in the same moving windows.
A point is declared an anomaly if its MAD value is more than k standard
deviations. The MAD is converted to a standard deviation using MAD * 1.4826,
which holds for normally distributed data.
The first bandwidth and last bandwidth observations cannot
be declared anomalies.
logical vector identifying which observations are anomalies.
Rob J Hyndman
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 10.2, https://OTexts.com/weird/.
set.seed(1) df <- tibble( time = seq(41), y = c(rnorm(20), 5, rnorm(20)) ) |> mutate(hampel = hampel_anomalies(y, bandwidth = 3, k = 4)) df |> ggplot(aes(x = time, y = y)) + geom_line() + geom_point(data = df |> filter(hampel), col = "red")set.seed(1) df <- tibble( time = seq(41), y = c(rnorm(20), 5, rnorm(20)) ) |> mutate(hampel = hampel_anomalies(y, bandwidth = 3, k = 4)) df |> ggplot(aes(x = time, y = y)) + geom_line() + geom_point(data = df |> filter(hampel), col = "red")
For a dist_kde object, determine which highest density region (HDR)
each observation falls in, for one or more coverage probabilities. Some
densities are multimodal, so the HDR for a given prob can consist of
several disjoint regions; the returned columns identify the specific
region (1, 2, ...) that each observation falls in, ordered from
lowest to highest along the first variable, with NA for observations
outside the HDR at that prob.
hdr_regions(object, prob)hdr_regions(object, prob)
object |
A |
prob |
A numeric vector of probabilities giving the HDR coverage (between 0 and 1). |
A tibble containing the data used to estimate object,
along with one additional integer
column per element of prob (named hdr_<100*prob>) showing which region
of the corresponding HDR each observation falls in.
Rob J Hyndman
dist_kde(oldfaithful$duration) |> hdr_regions(c(0.5, 0.95)) dist_kde(oldfaithful[, c("duration", "waiting")]) |> hdr_regions(0.90)dist_kde(oldfaithful$duration) |> hdr_regions(c(0.5, 0.95)) dist_kde(oldfaithful[, c("duration", "waiting")]) |> hdr_regions(0.90)
Compute a table of highest density regions (HDR) for a distributional object.
The HDRs are returned as a tibble with one row per interval and columns:
prob (giving the probability coverage),
density (the value of the density at the boundary of the HDR),
For one dimensional density functions, the tibble also has columns
lower (the lower ends of the intervals), and
upper (the upper ends of the intervals).
hdr_table(object, prob)hdr_table(object, prob)
object |
Distributional object such as that returned by |
prob |
Vector of probabilities giving the HDR coverage (between 0 and 1) |
A tibble
Rob J Hyndman
Hyndman, R J (1996) "Computing and Graphing Highest Density Regions", The American Statistician, 50(2), 120–126. https://robjhyndman.com/publications/hdr/
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 2.7, 3.4. https://OTexts.com/weird/.
# Univariate HDRs c(dist_normal(), dist_kde(c(rnorm(100), rnorm(100, 3, 1)))) |> hdr_table(c(0.5, 0.95)) dist_kde(oldfaithful$duration) |> hdr_table(0.95) # Bivariate HDRs dist_kde(oldfaithful[, c("duration", "waiting")]) |> hdr_table(0.90)# Univariate HDRs c(dist_normal(), dist_kde(c(rnorm(100), rnorm(100, 3, 1)))) |> hdr_table(c(0.5, 0.95)) dist_kde(oldfaithful$duration) |> hdr_table(0.95) # Bivariate HDRs dist_kde(oldfaithful[, c("duration", "waiting")]) |> hdr_table(0.90)
Bandwidth matrices are estimated in several ways including a normal reference rule, a robust version of the normal reference rule (default), a plugin estimator, a smoothed cross-validation estimator, or using the approach of Hyndman, Kandanaarachchi & Turner (2026). Details of each method are given in Hyndman (2026).
kde_bandwidth( data, method = c("robust", "normal", "plugin", "scv", "lookout"), ... )kde_bandwidth( data, method = c("robust", "normal", "plugin", "scv", "lookout"), ... )
data |
A numeric matrix or data frame. |
method |
A character string giving the method to use. Possibilities are:
|
... |
Additional arguments are ignored. |
A matrix of bandwidths (or a scalar in the case of univariate data).
Rob J Hyndman
Hyndman, R J, Kandanaarachchi, S & Turner, K (2026) "When lookout sees crackle: Anomaly detection via kernel density estimation", unpublished. https://robjhyndman.com/publications/lookout2.html
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 2.9 and 3.9, https://OTexts.com/weird/.
# Univariate bandwidth calculation kde_bandwidth(oldfaithful$duration) # Bivariate bandwidth calculation kde_bandwidth(oldfaithful[, c("duration", "waiting")])# Univariate bandwidth calculation kde_bandwidth(oldfaithful$duration) # Bivariate bandwidth calculation kde_bandwidth(oldfaithful[, c("duration", "waiting")])
Compute local outlier factors using k nearest neighbours. A local
outlier factor is a measure of how anomalous each observation is based on
the density of neighbouring points.
The function uses dbscan::lof to do the calculation.
lof_scores(y, k = 10, ...)lof_scores(y, k = 10, ...)
y |
Numerical matrix or vector of data |
k |
Number of neighbours to include, not counting the observation itself. Default: 10. |
... |
Additional arguments passed to |
Numerical vector containing LOF values. An observation has an infinite
LOF when its neighbourhood includes at least k + 1 identical observations
(whose local reachability density is infinite) but it is not one of them;
the identical observations themselves have LOF values of 1.
Rob J Hyndman
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 6.6, https://OTexts.com/weird/.
dbscan::lof
y <- c(rnorm(49), 5) lof_scores(y)y <- c(rnorm(49), 5) lof_scores(y)
A multivariate version of base::scale(), that takes account
of the covariance matrix of the data. By default, robust estimates are used:
the centers are removed using medians,
the scale function for univariate data is s_Qn,
and the covariance matrix for multivariate data is estimated using a robust MCD estimate.
The data are scaled using the Cholesky decomposition of
the inverse (co)variance. Then the scaled data are returned.
Details of the methods are provided by Hyndman (2026).
mvscale( object, center = stats::median, scale = robustbase::s_Qn, cov = robustbase::covMcd, alpha = 0.9, warning = TRUE, ... )mvscale( object, center = stats::median, scale = robustbase::s_Qn, cov = robustbase::covMcd, alpha = 0.9, warning = TRUE, ... )
object |
A vector, matrix, or data frame containing some numerical data. |
center |
A function to compute the center of each numerical variable. Set to NULL if no centering is required. |
scale |
A function to scale each numerical variable. When
|
cov |
A function to compute the covariance matrix. Set to NULL if no rotation required. |
alpha |
When |
warning |
Should a warning be issued if non-numeric columns are ignored? |
... |
Other arguments are passed to |
Optionally, the centering and scaling can be done for each variable
separately, by setting cov = NULL, so there is no rotation of the data,
Also optionally, non-robust methods can be used by specifying center = mean,
scale = stats::sd(), and cov = stats::cov(). Any non-numeric columns are retained
with a warning.
Missing values are removed before the centers, scale and cov are estimated.
A vector, matrix or data frame of the same size and class as object,
but with numerical variables replaced by scaled versions (renamed if they have been rotated).
Rob J Hyndman
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 2.6, 3.6 and 3.7.
base::scale(), stats::sd(), stats::cov(), robustbase::covMcd(), robustbase::covOGK(), robustbase::s_Qn()
# Univariate z-scores z <- mvscale(oldfaithful$duration, center = mean, scale = sd) # Non-robust scaling with no rotation oldfaithful |> mvscale(center = mean, scale = sd, cov = NULL, warning = FALSE) # Non-robust scaling with rotation oldfaithful |> mvscale(center = mean, scale = sd, cov = stats::cov, warning = FALSE) # Robust scaling and rotation oldfaithful |> mvscale(warning = FALSE)# Univariate z-scores z <- mvscale(oldfaithful$duration, center = mean, scale = sd) # Non-robust scaling with no rotation oldfaithful |> mvscale(center = mean, scale = sd, cov = NULL, warning = FALSE) # Non-robust scaling with rotation oldfaithful |> mvscale(center = mean, scale = sd, cov = stats::cov, warning = FALSE) # Robust scaling and rotation oldfaithful |> mvscale(warning = FALSE)
A synthetic data set containing 1000 observations on 10 variables generated from independent standard normal distributions.
n01n01
A data frame with 1000 rows and 10 columns.
Data frame
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 1.4, https://OTexts.com/weird/.
n01n01
A data set containing data on recorded eruptions of the Old Faithful Geyser in Yellowstone National Park, Wyoming, USA, from 14 January 2017 to 29 December 2023. Recordings are incomplete, especially during the winter months when observers may not be present.
oldfaithfuloldfaithful
A data frame with 2097 rows and 4 columns:
Time eruption started
Duration of eruption as recorded
Duration of eruption in seconds
Time to the following eruption in seconds
Data frame
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 1.4, https://OTexts.com/weird/.
oldfaithful |> ggplot(aes(x = duration, y = waiting)) + geom_point()oldfaithful |> ggplot(aes(x = duration, y = waiting)) + geom_point()
Draw an outlier map showing the score distance and orthogonal
distance of each observation from a projection or principal component
analysis. The score distance measures how far an observation lies from the
centre within the projection subspace, while the orthogonal distance
measures how far it lies from the subspace. Pass object (the output of
stats::prcomp() or an rrcov Pca* function); otherwise supply scores and
loadings together with the original data. For a prcomp object, use its
rank. argument to set the number of retained components.
When object is a PCA-like object and show_thresholds = TRUE, the
score-distance and orthogonal-distance cutoffs are drawn as dashed lines and
observations are coloured by type:
small score and orthogonal distance.
large score distance, small orthogonal distance.
small score distance, large orthogonal distance.
large score and orthogonal distance.
The cutoffs are only defined for PCA-like objects, so show_thresholds is
ignored when scores and loadings are passed directly.
outlier_map( object = NULL, data = NULL, scores = NULL, loadings = NULL, show_thresholds = TRUE, ... )outlier_map( object = NULL, data = NULL, scores = NULL, loadings = NULL, show_thresholds = TRUE, ... )
object |
Optionally, the output of |
data |
The original data matrix or data frame used to compute the
projection, scaled if the projection was computed on scaled data. This is
required to compute the orthogonal distances, except when |
scores |
A matrix or data frame of scores, with one column per retained
component. Ignored if |
loadings |
A matrix or data frame of loadings, with one column per
retained component. Ignored if |
show_thresholds |
If |
... |
Additional arguments passed to |
A ggplot object.
Rob J Hyndman
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 9, https://OTexts.com/weird/highdim.html.
Y <- oldfaithful[, c("duration", "waiting")] prcomp(Y, scale = TRUE, rank. = 1) |> outlier_map(data = Y)Y <- oldfaithful[, c("duration", "waiting")] prcomp(Y, scale = TRUE, rank. = 1) |> outlier_map(data = Y)
Peirce's criterion and Chauvenet's criterion were both proposed in the 1800s as a way of determining what observations should be rejected in a univariate sample.
peirce_anomalies(y) chauvenet_anomalies(y)peirce_anomalies(y) chauvenet_anomalies(y)
y |
numerical vector of observations |
These functions take a univariate sample y and return a logical
vector indicating which observations should be considered anomalies according
to either Peirce's criterion or Chauvenet's criterion.
A logical vector
Rob J Hyndman
Peirce, B (1852). Criterion for the rejection of doubtful observations. The Astronomical Journal, 2(21), 161–163.
Chauvenet, W (1863). 'Method of least squares'. Appendix to Manual of Spherical and Practical Astronomy, Vol.2, Lippincott, Philadelphia, pp.469-566.
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Section 4.3, https://OTexts.com/weird/.
y <- rnorm(1000) tibble(y = y) |> filter(peirce_anomalies(y)) tibble(y = y) |> filter(chauvenet_anomalies(y))y <- rnorm(1000) tibble(y = y) |> filter(peirce_anomalies(y)) tibble(y = y) |> filter(chauvenet_anomalies(y))
Every death from a road crash in Australia from January 1989 to December 2025,
from the Australian Road Deaths Database (ARDD). Each row is a person killed;
people killed in the same crash share a crash_id. A death is counted if the
crash occurred on a public road, was unintentional, and the person died within
30 days.
road_deathsroad_deaths
A data frame with 58168 rows and 20 columns:
Crash identifier
State or territory where the crash occurred
Year of crash
Month of crash (1–12)
Day of the week of the crash
Time of crash in hours after midnight (e.g. 14.5 is 2:30pm)
Single or multiple vehicle crash
Was a bus involved?
Was a heavy rigid truck involved?
Was an articulated truck involved?
Posted speed limit at the crash location (km/h)
Road user type of the person killed
Sex of the person killed
Age of the person killed (years)
ABS remoteness area of the crash location (ASGS 2021)
ABS Statistical Area Level 4 of the crash location (ASGS 2021)
Local government area of the crash location (ASGS 2021)
Type of road
Did the crash occur in the 12 days from 23 December?
Did the crash occur in the 5 days from the Thursday before Good Friday?
Missing values (coded as -9 or "Unknown" in the source) are NA. Some
variables were not recorded for early years: heavy_rigid_truck is mostly
missing before 2002, and remoteness, sa4, lga and road_type are
missing before 2014 and incomplete until 2017.
Data frame
Bureau of Infrastructure and Transport Research Economics (2026). Australian Road Deaths Database, fatalities. Downloaded 9 October 2026. Licensed under CC BY. https://catalogue.data.infrastructure.gov.au/dataset/australian-road-deaths-database
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 10, https://OTexts.com/weird/text.html.
road_deaths |> ggplot(aes(x = age, fill = road_user)) + geom_histogram(binwidth = 1)road_deaths |> ggplot(aes(x = age, fill = road_user)) + geom_histogram(binwidth = 1)
Test if observations are anomalies according to the stray algorithm.
stray_anomalies(y, ...)stray_anomalies(y, ...)
y |
A vector, matrix, or data frame consisting of numerical variables. |
... |
Other arguments are passed to |
Numerical vector containing logical values indicating if the observation is identified as an anomaly using the stray algorithm.
Rob J Hyndman
Talagala, P D, Hyndman, R J, and Smith-Miles, K (2021) Anomaly detection in high-dimensional data, Journal of Computational and Graphical Statistics, 30(2), 360-374.
# Univariate data y <- c(6, rnorm(49)) stray_anomalies(y) # Bivariate data y <- cbind(rnorm(50), c(5, rnorm(49))) stray_anomalies(y)# Univariate data y <- c(6, rnorm(49)) stray_anomalies(y) # Bivariate data y <- cbind(rnorm(50), c(5, rnorm(49))) stray_anomalies(y)
Compute stray scores indicating how anomalous each observation is.
stray_scores(y, ...)stray_scores(y, ...)
y |
A vector, matrix, or data frame consisting of numerical variables. |
... |
Other arguments are passed to |
Numerical vector containing stray scores.
Rob J Hyndman
Talagala, P D, Hyndman, R J, and Smith-Miles, K (2021) Anomaly detection in high-dimensional data, Journal of Computational and Graphical Statistics, 30(2), 360-374.
# Univariate data y <- c(6, rnorm(49)) scores <- stray_scores(y) # Indices of anomalies stray::find_threshold(scores, alpha = 0.01, outtail = "max", p = 0.5, tn = 50)# Univariate data y <- c(6, rnorm(49)) scores <- stray_scores(y) # Indices of anomalies stray::find_threshold(scores, alpha = 0.01, outtail = "max", p = 0.5, tn = 50)
A surprisal is given by where is the
density or probability mass function of the estimated or assumed distribution,
and is an observation. This is returned by surprisals().
A surprisal probability is the probability of a surprisal at least as extreme
as . This is returned by surprisals_prob()
surprisals(object, ...) surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, ... )surprisals(object, ...) surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, ... )
object |
A model or numerical data set |
... |
Other arguments are passed to the appropriate method. |
approximation |
Character string specifying the method to use in computing the surprisal probabilities. See Details below. |
threshold_probability |
Probability threshold when computing the GPD
approximation. This is the probability below which the GPD is fitted. Only
used if |
The surprisal probabilities may be computed in three different ways.
When approximation = "none" (the default), the surprisal probabilities are computed
using the same distribution that was used to compute the surprisal values.
Under this option, surprisal probabilities are equal to 1 minus the
coverage probability of the largest HDR that contains each value. Surprisal
probabilities smaller than 1e-6 are returned as 1e-6.
When approximation = "gdp", the surprisal probabilities are
computed using a Generalized Pareto Distribution fitted to the most extreme
surprisal values (those with probability less than threshold_probability).
For surprisal probabilities greater than threshold_probability, the value of
threshold_probability is returned. Under this option, the distribution is
used for computing the surprisal values but not for determining their
probabilities. Due to extreme value theory, the resulting probabilities should
be relatively insensitive to the distribution used in computing the surprisal
values.
When approximation = "empirical" (or "rank"), the surprisal probability of each
observation is estimated using the proportion of observations with
greater or equal surprisal values. This is a nonparametric approach that is also
insensitive to the distribution used in computing the surprisal values.
A numerical vector containing the surprisals or surprisal probabilities.
Rob J Hyndman
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 7, https://OTexts.com/weird/density.html.
Hyndman, R J & Frazier, D T (2026) "Anomaly detection using surprisals", https://robjhyndman.com/publications/surprisals.html.
For specific methods, see surprisals.numeric() and surprisals.lm(),
A surprisal is given by where is the
density or probability mass function of the estimated or assumed distribution,
and is an observation. This is returned by surprisals().
A surprisal probability is the probability of a surprisal at least as extreme
as . This is returned by surprisals_prob()
## S3 method for class 'lm' surprisals(object, loo = FALSE, ...) ## S3 method for class 'glm' surprisals(object, ...) ## S3 method for class 'gam' surprisals(object, ...) ## S3 method for class 'lm' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, loo = FALSE, ... ) ## S3 method for class 'glm' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, ... ) ## S3 method for class 'gam' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, ... )## S3 method for class 'lm' surprisals(object, loo = FALSE, ...) ## S3 method for class 'glm' surprisals(object, ...) ## S3 method for class 'gam' surprisals(object, ...) ## S3 method for class 'lm' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, loo = FALSE, ... ) ## S3 method for class 'glm' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, ... ) ## S3 method for class 'gam' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, ... )
object |
A model object such as returned by |
loo |
Should leave-one-out surprisals be computed? For computational
reasons, this is only available for |
... |
Other arguments are ignored. |
approximation |
Character string specifying the method to use in computing the surprisal probabilities. See Details below. |
threshold_probability |
Probability threshold when computing the GPD
approximation. This is the probability below which the GPD is fitted. Only
used if |
The surprisal probabilities may be computed in three different ways.
When approximation = "none" (the default), the surprisal probabilities are computed
using the same distribution that was used to compute the surprisal values.
Under this option, surprisal probabilities are equal to 1 minus the
coverage probability of the largest HDR that contains each value. Surprisal
probabilities smaller than 1e-6 are returned as 1e-6.
When approximation = "gdp", the surprisal probabilities are
computed using a Generalized Pareto Distribution fitted to the most extreme
surprisal values (those with probability less than threshold_probability).
For surprisal probabilities greater than threshold_probability, the value of
threshold_probability is returned. Under this option, the distribution is
used for computing the surprisal values but not for determining their
probabilities. Due to extreme value theory, the resulting probabilities should
be relatively insensitive to the distribution used in computing the surprisal
values.
When approximation = "empirical" (or "rank"), the surprisal probability of each
observation is estimated using the proportion of observations with
greater or equal surprisal values. This is a nonparametric approach that is also
insensitive to the distribution used in computing the surprisal values.
A numerical vector containing the surprisals or surprisal probabilities.
Rob J Hyndman
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 7, https://OTexts.com/weird/density.html.
Hyndman, R J & Frazier, D T (2026) "Anomaly detection using surprisals", https://robjhyndman.com/publications/surprisals.html.
For specific methods, see surprisals.numeric() and surprisals.lm(),
# A linear model (i.e., a conditional Gaussian distribution) lm_of <- lm(waiting ~ duration, data = oldfaithful) oldfaithful |> mutate( fscore = surprisals_prob(lm_of), prob = surprisals_prob(lm_of, loo = TRUE), ) |> ggplot(aes( x = duration, y = waiting, color = prob < 0.01 )) + geom_point() # A Poisson GLM glm_breaks <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson) warpbreaks |> mutate(prob = surprisals_prob(glm_breaks)) |> filter(prob < 0.05)# A linear model (i.e., a conditional Gaussian distribution) lm_of <- lm(waiting ~ duration, data = oldfaithful) oldfaithful |> mutate( fscore = surprisals_prob(lm_of), prob = surprisals_prob(lm_of, loo = TRUE), ) |> ggplot(aes( x = duration, y = waiting, color = prob < 0.01 )) + geom_point() # A Poisson GLM glm_breaks <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson) warpbreaks |> mutate(prob = surprisals_prob(glm_breaks)) |> filter(prob < 0.05)
A surprisal is given by where is the
density or probability mass function of the estimated or assumed distribution,
and is an observation. This is returned by surprisals().
A surprisal probability is the probability of a surprisal at least as extreme
as . This is returned by surprisals_prob()
## S3 method for class 'numeric' surprisals(object, distribution = dist_kde(object, ...), loo = FALSE, ...) ## S3 method for class 'matrix' surprisals(object, distribution = dist_kde(object, ...), loo = FALSE, ...) ## S3 method for class 'data.frame' surprisals(object, distribution = dist_kde(object, ...), loo = FALSE, ...) ## S3 method for class 'numeric' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, distribution = dist_kde(object, ...), loo = FALSE, ... ) ## S3 method for class 'matrix' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, distribution = dist_kde(object, ...), loo = FALSE, ... ) ## S3 method for class 'data.frame' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, distribution = dist_kde(object, ...), loo = FALSE, ... )## S3 method for class 'numeric' surprisals(object, distribution = dist_kde(object, ...), loo = FALSE, ...) ## S3 method for class 'matrix' surprisals(object, distribution = dist_kde(object, ...), loo = FALSE, ...) ## S3 method for class 'data.frame' surprisals(object, distribution = dist_kde(object, ...), loo = FALSE, ...) ## S3 method for class 'numeric' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, distribution = dist_kde(object, ...), loo = FALSE, ... ) ## S3 method for class 'matrix' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, distribution = dist_kde(object, ...), loo = FALSE, ... ) ## S3 method for class 'data.frame' surprisals_prob( object, approximation = c("none", "gpd", "empirical", "rank"), threshold_probability = 0.1, distribution = dist_kde(object, ...), loo = FALSE, ... )
object |
A numerical data set (either a vector, matrix, or a data.frame containing only numerical columns). |
distribution |
A distribution object. By default, a kernel density
estimate is computed from the data |
loo |
Should leave-one-out surprisals be computed? |
... |
Other arguments are passed to the appropriate method. |
approximation |
Character string specifying the method to use in computing the surprisal probabilities. See Details below. For a multivariate data set, it needs to be set to either "gpd" or "empirical". |
threshold_probability |
Probability threshold when computing the GPD
approximation. This is the probability below which the GPD is fitted. Only
used if |
The surprisal probabilities may be computed in three different ways.
When approximation = "none" (the default), the surprisal probabilities are computed
using the same distribution that was used to compute the surprisal values.
Under this option, surprisal probabilities are equal to 1 minus the
coverage probability of the largest HDR that contains each value. Surprisal
probabilities smaller than 1e-6 are returned as 1e-6.
When approximation = "gdp", the surprisal probabilities are
computed using a Generalized Pareto Distribution fitted to the most extreme
surprisal values (those with probability less than threshold_probability).
For surprisal probabilities greater than threshold_probability, the value of
threshold_probability is returned. Under this option, the distribution is
used for computing the surprisal values but not for determining their
probabilities. Due to extreme value theory, the resulting probabilities should
be relatively insensitive to the distribution used in computing the surprisal
values.
When approximation = "empirical" (or "rank"), the surprisal probability of each
observation is estimated using the proportion of observations with
greater or equal surprisal values. This is a nonparametric approach that is also
insensitive to the distribution used in computing the surprisal values.
A numerical vector containing the surprisals or surprisal probabilities.
Rob J Hyndman
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 7, https://OTexts.com/weird/density.html.
Hyndman, R J & Frazier, D T (2026) "Anomaly detection using surprisals", https://robjhyndman.com/publications/surprisals.html.
# Univariate data tibble( y = c(5, rnorm(49)), p_kde = surprisals_prob(y, loo = TRUE), p_normal = surprisals_prob(y, distribution = dist_normal()), p_zscore = 2 * (1 - pnorm(abs(y))) ) tibble( y = n01$v1, prob1 = surprisals_prob(y), prob2 = surprisals_prob(y, loo = TRUE), prob3 = surprisals_prob(y, distribution = dist_normal()), prob4 = surprisals_prob(y, distribution = dist_normal(), approximation = "gpd") ) |> arrange(prob1) # Bivariate data tibble( x = rnorm(50), y = c(5, rnorm(49)), prob = surprisals_prob(cbind(x, y), approximation = "gpd") ) oldfaithful |> mutate( s = surprisals(cbind(duration, waiting), loo = TRUE), p = surprisals_prob(cbind(duration, waiting), loo = TRUE, approximation = "gpd") ) |> arrange(p)# Univariate data tibble( y = c(5, rnorm(49)), p_kde = surprisals_prob(y, loo = TRUE), p_normal = surprisals_prob(y, distribution = dist_normal()), p_zscore = 2 * (1 - pnorm(abs(y))) ) tibble( y = n01$v1, prob1 = surprisals_prob(y), prob2 = surprisals_prob(y, loo = TRUE), prob3 = surprisals_prob(y, distribution = dist_normal()), prob4 = surprisals_prob(y, distribution = dist_normal(), approximation = "gpd") ) |> arrange(prob1) # Bivariate data tibble( x = rnorm(50), y = c(5, rnorm(49)), prob = surprisals_prob(cbind(x, y), approximation = "gpd") ) oldfaithful |> mutate( s = surprisals(cbind(duration, waiting), loo = TRUE), p = surprisals_prob(cbind(duration, waiting), loo = TRUE, approximation = "gpd") ) |> arrange(p)
Weekly deaths and death rates in the USA, by sex and age group, from the second week of 2015 to the last week of 2025. The data are from the Short-Term Mortality Fluctuations (STMF) series of the Human Mortality Database.
us_mortalityus_mortality
A data frame with 10314 rows and 6 columns:
Year
ISO week of the year
"Female", "Male" or "Total"
Age group: "0-14", "15-64", "65-74", "75-84", "85+" or "Total"
Number of deaths
Weekly death rate: deaths divided by population exposure
Data frame
Human Mortality Database. Max Planck Institute for Demographic Research (Germany), University of California, Berkeley (USA), and French Institute for Demographic Studies (France). https://www.mortality.org
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 11, https://OTexts.com/weird/timeseries.html.
us_mortality |> filter(Sex == "Total", Age != "Total") |> mutate(time = Year + (Week - 1) / 52) |> ggplot(aes(x = time, y = Mortality, colour = Age)) + geom_line() + scale_y_log10()us_mortality |> filter(Sex == "Total", Age != "Total") |> mutate(time = Year + (Week - 1) / 52) |> ggplot(aes(x = time, y = Mortality, colour = Age)) + geom_line() + scale_y_log10()
Co-voting networks for the US Senate from the 40th to the 113th Congress (1867 to 2015). Each Congress gives an undirected network in which each node is a senator, and an edge connects two senators if they voted the same way (both yea or both nay) on at least 80% of the bills for which they were both present. The data are taken from Lee, Li and Wilson (2020). Although they state that Independent senators were excluded, a few Independent and minor-party senators are included, coded as Democrat or Republican.
us_senate_edges us_senate_membersus_senate_edges us_senate_members
us_senate_edges is a data frame with 73802 rows and
3 columns:
Congress number
Node number of one senator
Node number of the other senator
us_senate_members is a data frame with 7272 rows and
6 columns:
Congress number
Node number of the senator
Name of the senator (surname first)
Two-letter state abbreviation ("USA" for the President)
"Democrat" or "Republican"
ICPSR identifier of the senator, as used by Voteview
An object of class tbl_df (inherits from tbl, data.frame) with 7272 rows and 6 columns.
us_senate_edges contains the edges of every network, and us_senate_members
contains the name, state and party of each senator. Senators are numbered
from 1 within each Congress, so the same number refers to different
senators in different Congresses; use icpsr to follow a senator across
Congresses. From the 84th Congress onwards, node 1 is the President, whose
announced positions on bills are treated as votes.
Data frame
Lee, J., Li, G., & Wilson, J. D. (2020). Varying-coefficient models for dynamic networks. Computational Statistics & Data Analysis, 152, 107052. doi:10.1016/j.csda.2020.107052. Data available from https://github.com/jihuilee/VCERGM.
Names, states and ICPSR identifiers from Lewis, J. B., Poole, K., Rosenthal, H., Boche, A., Rudkin, A., & Sonnet, L. Voteview: Congressional roll-call votes database. https://voteview.com/
Hyndman, R J (2026) "That's weird: Anomaly detection using R", Chapter 15, https://OTexts.com/weird/network.html.
# Number of co-voting edges in each Congress us_senate_edges |> count(congress) # Party composition of the 100th Congress us_senate_members |> filter(congress == 100) |> count(party) # Congresses in which Patrick Leahy served us_senate_members |> filter(icpsr == 14307)# Number of co-voting edges in each Congress us_senate_edges |> count(congress) # Party composition of the 100th Congress us_senate_members |> filter(congress == 100) |> count(party) # Congresses in which Patrick Leahy served us_senate_members |> filter(icpsr == 14307)