Three types of surprise

Weird times: unusual observations in a historical series

Weird series: an unusual series in a collection

Surveillance: an unusual observation just arrived

Outline

  1. Surprisals
  2. Weird times
  3. Weird series
  4. Surveillance
  5. Summary

Surprisals

01

  1. Surprisals
  2. Weird times
  3. Weird series
  4. Surveillance
  5. Summary

Surprisal probabilities

The probability of an observation at least as unlikely as y_t under density f:

p_t = \mathbb{P}\big(f(Y) \le f(y_t)\big), \qquad Y \sim f

An observation is an anomaly with respect to f if p_t < \alpha.

Surprisals

The surprisal of an observation y_t under density f is s_t = -\log f(y_t)

Surprisals are already familiar:

Forecasting

Log score for evaluating density forecasts

Information theory

Expected surprisal = entropy

Likelihood

Sum of surprisals = -log likelihood

Let S = -\log f(Y) where Y \sim f, and G(s) = \mathbb{P}(S \le s).
Then p_t = \mathbb{P}(S \ge s_t) = 1 - G(s_t).

An anomaly is an extreme value of the surprisal distribution.

fmay be complicated

Multimodal, multivariate or conditional: G is still a univariate distribution on \mathbb{R}.

fmay be wrong

The upper tail probabilities of G can be estimated using EVT.

Surprisals have GPD tails

Theorem: surprisals have Generalised Pareto Distribution (GPD) tails

If the upper tail of S = -\log f(Y) is smooth, then as the threshold u grows,

\mathbb{P}\left( \frac{S - u}{a(u)} \le y \;\middle|\; S > u \right) \longrightarrow \begin{cases} 1 - (1 + \xi y)_+^{-1/\xi}, & \xi \ne 0,\\ 1 - e^{-y}, & \xi = 0. \end{cases}

\xi < 0: finite endpoint

S bounded above
e.g., truncated Normal data

\xi = 0: exponential tail

\mathbb{E}\,e^{\lambda S} < \infty
e.g., Normal or t data

\xi > 0: Pareto tail

Only \mathbb{E}\,S^q < \infty, with \xi \le 1/q
e.g., t data, Normal f

Conditions are on S, not Y: they hold even if f is wrong.

Approximate GPD surprisal probabilities

Only needs f to get the ordering of the extreme surprisals correct.

Univariate experiment

What happens when f is mis-specified?

Data

1000 observations from N(0,1), and 1000 from t(4)

Surprisals

Computed under N(0,1) and under t(4)

Probabilities

From f itself, or from a GPD fitted to the surprisals

Univariate experiment

Weird times

02

  1. Surprisals
  2. Weird times
  3. Weird series
  4. Surveillance
  5. Summary

French mortality rates

library(weird)
fr_mortality
# A tibble: 31,648 × 4
   Year   Age Sex    Mortality
  <int> <int> <chr>      <dbl>
1  1816     0 Female    0.187 
2  1816     1 Female    0.0467
3  1816     2 Female    0.0339
4  1816     3 Female    0.0229
5  1816     4 Female    0.0160
6  1816     5 Female    0.0138
7  1816     6 Female    0.0121
8  1816     7 Female    0.0104
# ℹ 31,640 more rows
  • Annual data: 1816–1999
  • Ages 0–85, by sex

French mortality rates

Hampel identifier

m_t = \text{median}(y_{t-w}, \dots, y_{t+w})

a_t = \text{median}(|y_{t-w} - m_t|, \dots, |y_{t+w} - m_t|)

  • Anomaly if |y_t - m_t| > ka_t
  • Implicit model: y_t \sim \text{N}(m_t, \sigma_t^2) with \sigma_t = a_t/0.6745
  • Better: compute surprisals under Normal, with probabilities from GPD

French mortality rates

fr_anomalies <- fr_mortality |>
  group_by(Age, Sex) |>
  mutate(
    hampel = hampel_anomalies(
      Mortality,
      bandwidth = 7,
      alpha = 0.01,
      approximation = "gpd"
    )
  ) |>
  ungroup() |>
  filter(hampel)
fr_anomalies |>
  count(Year, Sex) |>
  filter(n >= 3)
# A tibble: 328 × 5
    Year   Age Sex    Mortality hampel
   <int> <int> <chr>      <dbl> <lgl> 
 1  1832    40 Female    0.0167 TRUE  
 2  1832    42 Female    0.0171 TRUE  
 3  1832    45 Female    0.0174 TRUE  
 4  1832    47 Female    0.0191 TRUE  
 5  1832    49 Female    0.0206 TRUE  
 6  1832    50 Female    0.0217 TRUE  
 7  1832    52 Female    0.0242 TRUE  
 8  1832    53 Female    0.0253 TRUE  
 9  1832    54 Female    0.0264 TRUE  
10  1832    55 Female    0.0277 TRUE  
11  1832    59 Female    0.0365 TRUE  
12  1832    60 Female    0.0403 TRUE  
13  1832    61 Female    0.0445 TRUE  
14  1832    62 Female    0.0481 TRUE  
15  1832    64 Female    0.0531 TRUE  
# ℹ 313 more rows
# A tibble: 18 × 3
    Year Sex        n
   <int> <chr>  <int>
 1  1832 Female    21
 2  1832 Male      14
 3  1834 Female     4
 4  1834 Male       5
 5  1837 Female     9
 6  1837 Male       5
 7  1849 Female    14
 8  1849 Male       4
 9  1870 Female     4
10  1870 Male       8
11  1871 Female    64
12  1871 Male      51
13  1914 Male       8
14  1915 Male      12
15  1918 Female    35
16  1918 Male       9
17  1940 Female    11
18  1940 Male      11

French mortality rates

US weekly mortality rates

us_stmf
# A tsibble: 2,865 x 3 [1W]
# Key:       Age group [5]
      Week `Age group` Mortality
    <week> <chr>           <dbl>
1 2015 W02 0-14         0.000522
2 2015 W03 0-14         0.000537
3 2015 W04 0-14         0.000532
4 2015 W05 0-14         0.000526
5 2015 W06 0-14         0.000494
6 2015 W07 0-14         0.000528
7 2015 W08 0-14         0.000530
8 2015 W09 0-14         0.000522
# ℹ 2,857 more rows
  • Weekly data: 2015 W02 to 2025 W52
  • 5 age groups: 0-14, 15-64, 65-74, 75-84, 85+

US weekly mortality rates

Surprisals from an STL decomposition

Local medians fail with seasonality.

Use a better signal estimate:

  • STL: y_t = T_t + S_t + R_t (Cleveland et al., 1990; Bandara et al., 2025)
  • a: robust scale of R_t in each age group
  • e_t = R_t / a; surprisals under N(0,1), probabilities from GPD
stmf_scores <- us_stmf |>
  model(stl = STL(Mortality)) |> augment() |>
  group_by(`Age group`) |>
  mutate(
    a = robustbase::s_Qn(.resid),
    e = .resid / a
  ) |>
  ungroup() |>
  mutate(prob = surprisals_prob(e, "gpd", distribution = dist_normal())) |>
  arrange(prob)
# A tsibble: 2,865 x 4 [1W]
# Key:       Age group [5]
       Week `Age group` Mortality  .fitted
     <week> <chr>           <dbl>    <dbl>
 1 2015 W02 0-14         0.000522 0.000536
 2 2015 W02 15-64        0.00350  0.00352 
 3 2015 W02 65-74        0.0200   0.0199  
 4 2015 W02 75-84        0.0552   0.0534  
 5 2015 W02 85+          0.190    0.178   
 6 2015 W03 0-14         0.000537 0.000527
 7 2015 W03 15-64        0.00347  0.00349 
 8 2015 W03 65-74        0.0197   0.0200  
 9 2015 W03 75-84        0.0545   0.0535  
10 2015 W03 85+          0.189    0.181   
# ℹ 2,855 more rows

# A tsibble: 2,865 x 5 [1W]
# Key:       Age group [5]
       Week `Age group` Mortality  .fitted         a
     <week> <chr>           <dbl>    <dbl>     <dbl>
 1 2015 W02 0-14         0.000522 0.000536 0.0000202
 2 2015 W02 15-64        0.00350  0.00352  0.000123 
 3 2015 W02 65-74        0.0200   0.0199   0.000649 
 4 2015 W02 75-84        0.0552   0.0534   0.00156  
 5 2015 W02 85+          0.190    0.178    0.00518  
 6 2015 W03 0-14         0.000537 0.000527 0.0000202
 7 2015 W03 15-64        0.00347  0.00349  0.000123 
 8 2015 W03 65-74        0.0197   0.0200   0.000649 
 9 2015 W03 75-84        0.0545   0.0535   0.00156  
10 2015 W03 85+          0.189    0.181    0.00518  
# ℹ 2,855 more rows

# A tsibble: 2,865 x 6 [1W]
# Key:       Age group [5]
       Week `Age group` Mortality  .fitted         a      e
     <week> <chr>           <dbl>    <dbl>     <dbl>  <dbl>
 1 2015 W02 0-14         0.000522 0.000536 0.0000202 -0.684
 2 2015 W02 15-64        0.00350  0.00352  0.000123  -0.117
 3 2015 W02 65-74        0.0200   0.0199   0.000649   0.127
 4 2015 W02 75-84        0.0552   0.0534   0.00156    1.14 
 5 2015 W02 85+          0.190    0.178    0.00518    2.25 
 6 2015 W03 0-14         0.000537 0.000527 0.0000202  0.493
 7 2015 W03 15-64        0.00347  0.00349  0.000123  -0.163
 8 2015 W03 65-74        0.0197   0.0200   0.000649  -0.529
 9 2015 W03 75-84        0.0545   0.0535   0.00156    0.601
10 2015 W03 85+          0.189    0.181    0.00518    1.58 
# ℹ 2,855 more rows

# A tsibble: 2,865 x 7 [1W]
# Key:       Age group [5]
       Week `Age group` Mortality  .fitted         a      e  prob
     <week> <chr>           <dbl>    <dbl>     <dbl>  <dbl> <dbl>
 1 2015 W02 0-14         0.000522 0.000536 0.0000202 -0.684   0.1
 2 2015 W02 15-64        0.00350  0.00352  0.000123  -0.117   0.1
 3 2015 W02 65-74        0.0200   0.0199   0.000649   0.127   0.1
 4 2015 W02 75-84        0.0552   0.0534   0.00156    1.14    0.1
 5 2015 W02 85+          0.190    0.178    0.00518    2.25    0.1
 6 2015 W03 0-14         0.000537 0.000527 0.0000202  0.493   0.1
 7 2015 W03 15-64        0.00347  0.00349  0.000123  -0.163   0.1
 8 2015 W03 65-74        0.0197   0.0200   0.000649  -0.529   0.1
 9 2015 W03 75-84        0.0545   0.0535   0.00156    0.601   0.1
10 2015 W03 85+          0.189    0.181    0.00518    1.58    0.1
# ℹ 2,855 more rows

# A tsibble: 2,865 x 7 [1W]
# Key:       Age group [5]
       Week `Age group` Mortality .fitted        a     e    prob
     <week> <chr>           <dbl>   <dbl>    <dbl> <dbl>   <dbl>
 1 2020 W16 85+           0.202   0.154   0.00518   9.22 0.00174
 2 2020 W15 75-84         0.0617  0.0480  0.00156   8.77 0.00208
 3 2020 W15 85+           0.200   0.155   0.00518   8.71 0.00213
 4 2021 W36 15-64         0.00543 0.00441 0.000123  8.30 0.00251
 5 2020 W15 65-74         0.0248  0.0194  0.000649  8.23 0.00259
 6 2020 W16 75-84         0.0600  0.0472  0.00156   8.20 0.00262
 7 2020 W17 85+           0.193   0.151   0.00518   8.00 0.00285
 8 2021 W35 15-64         0.00538 0.00441 0.000123  7.89 0.00299
 9 2021 W37 15-64         0.00528 0.00435 0.000123  7.60 0.00341
10 2020 W53 75-84         0.0696  0.0581  0.00156   7.35 0.00382
# ℹ 2,855 more rows

US weekly mortality rates

US weekly mortality rates

Weird series

03

  1. Surprisals
  2. Weird times
  3. Weird series
  4. Surveillance
  5. Summary

Australian Pharmaceutical Benefits Scheme

PBS <- tsibbledata::PBS |>
  update_tsibble(key = c(Concession, Type, ATC2)) |>
  select(Month, Concession, Type, ATC2, Scripts)
PBS
# A tsibble: 67,596 x 5 [1M]
# Key:       Concession, Type, ATC2 [336]
      Month Concession   Type        ATC2  Scripts
      <mth> <chr>        <chr>       <chr>   <dbl>
 1 1991 Jul Concessional Co-payments A01     18228
 2 1991 Aug Concessional Co-payments A01     15327
 3 1991 Sep Concessional Co-payments A01     14775
 4 1991 Oct Concessional Co-payments A01     15380
 5 1991 Nov Concessional Co-payments A01     14371
 6 1991 Dec Concessional Co-payments A01     15028
 7 1992 Jan Concessional Co-payments A01     11040
 8 1992 Feb Concessional Co-payments A01     15165
 9 1992 Mar Concessional Co-payments A01     16898
10 1992 Apr Concessional Co-payments A01     18141
# ℹ 67,586 more rows

How do we compare 336 time series?

Feature: strength of trend

Feature: strength of seasonality

Feature space

A weird series is an anomaly in feature space. Features are mostly scale-free, so series are flagged for their shape and dynamics, not their size.

Australian Pharmaceutical Benefits Scheme

PBS_feat <- PBS |>
  features(Scripts, feature_set(pkgs = "feasts"))
# A tibble: 336 × 51
  Concession   Type        ATC2  trend_strength seasonal_strength_year
  <chr>        <chr>       <chr>          <dbl>                  <dbl>
1 Concessional Co-payments A01            0.845                  0.918
2 Concessional Co-payments A02            0.970                  0.943
3 Concessional Co-payments A03            0.956                  0.931
# ℹ 333 more rows
# ℹ 46 more variables: seasonal_peak_year <dbl>, seasonal_trough_year <dbl>,
#   spikiness <dbl>, linearity <dbl>, curvature <dbl>, stl_e_acf1 <dbl>,
#   stl_e_acf10 <dbl>, acf1 <dbl>, acf10 <dbl>, diff1_acf1 <dbl>,
#   diff1_acf10 <dbl>, diff2_acf1 <dbl>, diff2_acf10 <dbl>, season_acf1 <dbl>,
#   pacf5 <dbl>, diff1_pacf5 <dbl>, diff2_pacf5 <dbl>, season_pacf <dbl>,
#   zero_run_mean <dbl>, nonzero_squared_cv <dbl>, zero_start_prop <dbl>, …

336 series, 48 features each; 3 series have missing features.

Australian Pharmaceutical Benefits Scheme

Australian Pharmaceutical Benefits Scheme

PBS_feat <- PBS_feat |> na.omit()
PBS_prcomp <- PBS_feat |>
  select(where(is.numeric)) |>
  prcomp(scale = TRUE, rank = 2) |>
  augment(PBS_feat)
PBS_prcomp |>
  select(Concession:ATC2, starts_with(".fitted"))
# A tibble: 333 × 5
  Concession   Type        ATC2  .fittedPC1 .fittedPC2
  <chr>        <chr>       <chr>      <dbl>      <dbl>
1 Concessional Co-payments A01         1.06     0.210 
2 Concessional Co-payments A02        -3.20     1.77  
3 Concessional Co-payments A03        -1.19     0.931 
4 Concessional Co-payments A04         1.21     0.0575
5 Concessional Co-payments A05        -1.44    -1.65  
# ℹ 328 more rows

Australian Pharmaceutical Benefits Scheme

J06: Concessional, Safety net
C05: General, Co-payments
S02: General, Co-payments
A15: General, Co-payments
C09: Concessional, Co-payments

Kernel density surprisals

KDE

\boldsymbol{y}_{i} \in \mathbb{R}^d, \quad \hat{f}(\boldsymbol{y}) = \frac{1}{n} \sum_{i=1}^n |\boldsymbol{H}|^{-1/2} K\big(\boldsymbol{H}^{-1/2}(\boldsymbol{y} - \boldsymbol{y}_i)\big)

  • K: spherically symmetric kernel; \boldsymbol{H}: d \times d bandwidth matrix.
  • Each observation inflates its own density estimate: use leave-one-out.

LOO KDE values

\hat{f}_{-i}(\boldsymbol{y}_i) = \frac{1}{n - 1} \big( n \hat{f}(\boldsymbol{y}_i) - |\boldsymbol{H}|^{-1/2} K(\boldsymbol{0}) \big)

  • Surprisals s_i = -\log \hat{f}_{-i}(\boldsymbol{y}_i); probabilities from the GPD.

Australian Pharmaceutical Benefits Scheme

J06: Concessional, Safety net
C05: General, Co-payments
S02: General, Co-payments
A15: General, Co-payments
C09: Concessional, Co-payments

Australian Pharmaceutical Benefits Scheme

PBS_prcomp |>
  mutate(
    prob = surprisals_prob(
      PBS_prcomp |> select(.fittedPC1, .fittedPC2),
      loo = TRUE,
      approximation = "gpd",
      method = "lookout"
    )
  ) |>
  select(Concession, Type, ATC2, .fittedPC1, .fittedPC2, prob) |>
  arrange(prob) |>
  print(n = 5)
# A tibble: 333 × 6
  Concession   Type        ATC2  .fittedPC1 .fittedPC2     prob
  <chr>        <chr>       <chr>      <dbl>      <dbl>    <dbl>
1 Concessional Safety net  J06         7.05     -17.9  0       
2 General      Co-payments C05         8.51     -23.0  0       
3 General      Co-payments S02         5.54     -11.5  0.000444
4 General      Co-payments A15        -1.15      -8.86 0.0200  
5 General      Co-payments P03         5.10      -9.70 0.0235  
# ℹ 328 more rows

Surveillance

04

  1. Surprisals
  2. Weird times
  3. Weird series
  4. Surveillance
  5. Summary

PBS scripts by ATC group

pbs <- PBS |>
  group_by(ATC2) |>
  summarise(Scripts = sum(Scripts) / 1e3) |>
  mutate(t = as.numeric(Month - min(Month) + 1))
# A tsibble: 17,016 x 4 [1M]
# Key:       ATC2 [84]
   ATC2     Month Scripts     t
   <chr>    <mth>   <dbl> <dbl>
 1 A01   1991 Jul    22.6     1
 2 A01   1991 Aug    20.4     2
 3 A01   1991 Sep    21.4     3
 4 A01   1991 Oct    23.7     4
 5 A01   1991 Nov    23.5     5
 6 A01   1991 Dec    26.3     6
 7 A01   1992 Jan    22.0     7
 8 A01   1992 Feb    16.4     8
 9 A01   1992 Mar    17.2     9
10 A01   1992 Apr    18.8    10
# ℹ 17,006 more rows

Surveillance using forecast distributions

ETS refitted monthly for every series. Points: p < 0.01, p < 0.001.

Surprisal probabilities as p-values

  • If f_{i,t+1|t} is the true forecast density, and the surprisal distribution is continuous, then p_{i,t+1} \mid y_{i,1},\dots,y_{i,t} \sim \text{U}(0,1) so the p_{i,t} are iid U(0,1) over t.
  • Thresholding at \alpha gives a false alarm rate of \alpha per series per period.
  • Model updated each period, so a level shift is flagged when it happens, then becomes the “new normal”.
  • ETS: Gaussian forecast distributions, so p_{i,t+1} is a two-sided z-test on the one-step forecast error.

A first iteration: t = 36

pbs_fc <- pbs |>
  filter(t <= 36) |>
  group_by_key() |>
  filter(cumsum(Scripts) > 0, n() >= 36) |>
  ungroup() |>
  model(ets = ETS(Scripts)) |>
  forecast(h = 1)
# A fable: 82 x 5 [1M]
# Key:     ATC2, .model [82]
  ATC2  .model    Month       Scripts  .mean
  <chr> <chr>     <mth>        <dist>  <dbl>
1 A01   ets    1994 Jul    N(23, 2.1)  23.1 
2 A02   ets    1994 Jul  N(590, 1054) 590.  
3 A03   ets    1994 Jul     N(84, 19)  84.2 
4 A04   ets    1994 Jul     N(69, 15)  69.4 
5 A06   ets    1994 Jul    N(33, 4.2)  33.1 
6 A07   ets    1994 Jul     N(74, 17)  73.6 
7 A09   ets    1994 Jul N(3.7, 0.029)   3.68
# ℹ 75 more rows

A first iteration: t = 36

pbs_fc |>
  rename(fcast = Scripts) |>
  left_join(pbs, by = c("ATC2", "Month")) |>
  mutate(prob = surprisals_prob(Scripts, distribution = fcast))
# A tsibble: 82 x 8 [?]
# Key:       ATC2, .model [82]
   ATC2  .model    Month         fcast  .mean Scripts     t   prob
   <chr> <chr>     <mth>        <dist>  <dbl>   <dbl> <dbl>  <dbl>
 1 A01   ets    1994 Jul    N(23, 2.1)  23.1    20.9     37 0.126 
 2 A02   ets    1994 Jul  N(590, 1054) 590.    516.      37 0.0233
 3 A03   ets    1994 Jul     N(84, 19)  84.2    80.5     37 0.393 
 4 A04   ets    1994 Jul     N(69, 15)  69.4    66.1     37 0.396 
 5 A06   ets    1994 Jul    N(33, 4.2)  33.1    29.2     37 0.0603
 6 A07   ets    1994 Jul     N(74, 17)  73.6    68.5     37 0.223 
 7 A09   ets    1994 Jul N(3.7, 0.029)   3.68    3.32    37 0.0313
 8 A10   ets    1994 Jul    N(166, 54) 166.    150.      37 0.0276
 9 A11   ets    1994 Jul      N(30, 3)  30.2    27.8     37 0.156 
10 A12   ets    1994 Jul    N(100, 23) 100.0    87.9     37 0.0121
# ℹ 72 more rows

Rolling forward through all series

267 observations with p_{i,t} < 0.001, out of 13983.

Surveillance anomalies

Surveillance across several series

Some events affect many series at once: mildly unusual in each, very unusual together.

  • M_t: number of series monitored at time t
  • N_t: number of series with p_{i,t} < 0.001
  • Independent series: N_t \sim \text{Binomial}(M_t, \pi)
  • Surprisal probability: \mathbb{P}(N_t \ge n_t)
  • \pi = 0.001 if forecast densities correct
  • ETS forecast densities too narrow: actual rate is 0.019
  • Solution: estimate \pi from all previous months

Surveillance across several series

Summary

05

  1. Surprisals
  2. Weird times
  3. Weird series
  4. Surveillance
  5. Summary

One idea, three types of surprise

Weird times
f: density of residuals from the signal

Weird series
f: kernel density of the features

Surveillance
f: one-step forecast density

  • Surprisal probabilities: p_t = \mathbb{P}\big(f(Y) \le f(y_t)\big)
  • GPD tail keeps p_t approximately right even when f is wrong

Coauthors

Sevvandi
Kandanaarachchi
CSIRO

Kate
Turner
ANU

David
Frazier
Monash U

Find out more

R package

weird on CRAN

install.packages("weird")

Book

That’s weird: Anomaly detection using R

In preparation at
OTexts.com/weird

Chapter 11: Time series data

Slides, links and related papers:
robjhyndman.com/uconn2026