Tutorial for sequential spatial ordinal pattern tests

Libraries

using StatsOrdinalPatterns, Random, Statistics, DataFrames, CairoMakie

Sequential testing with spatial ordinal patterns

Here, the monitored stream is a sequence of images: at each time \(t\) a new \(M \times N\) image arrives, all its \(2 \times 2\) windows at delays \((d_1, d_2)\) are classified into SOP types, and the resulting type frequencies \(\hat{\mathbf{p}}_t\) are smoothed exponentially (EWMA),

\[ \hat{\mathbf{p}}_t^{\text{EWMA}} = \lambda\,\hat{\mathbf{p}}_t + (1 - \lambda)\,\hat{\mathbf{p}}_{t-1}^{\text{EWMA}}, \qquad \hat{\mathbf{p}}_0^{\text{EWMA}} = (1/3, 1/3, 1/3). \]

The chart statistic is computed from the smoothed frequencies and the chart signals as soon as \(\lvert\text{stat}\rvert > cl\). Unlike the time series case, each single time point already contributes \((M-d_1)(N-d_2)\) patterns, so the chart reacts quickly.

Designing the chart means fixing \(\lambda\) and choosing cl such that the in-control ARL matches a target L0. See Adämmer et al. (2026).

Key functions

Alarm rules

The four Tau/Kappa statistics are used for monitoring; all of them are two-sided, so the chart signals when \(\lvert\text{stat}\rvert > cl\).

Statistic Constructor Definition Under \(H_0\)
\(\hat{\tau}\) TauHat() \(\hat{p}_1 - 1/3\) \(0\)
\(\hat{\kappa}\) KappaHat() \(\hat{p}_2 - \hat{p}_3\) \(0\)
\(\tilde{\tau}\) TauTilde() \(\tilde{p}_3 - 1/3\) \(0\)
\(\tilde{\kappa}\) KappaTilde() \(\tilde{p}_1 - \tilde{p}_2\) \(0\)

Setting refinement to RotationType(), DirectionType(), or DiagonalType() replaces the three classical SOP types by the corresponding six refined types. The refinement argument must be identical in cl_sop, arl_sop_ic, and arl_sop_oc — otherwise the critical limit does not belong to the monitored statistic.

Setup

M, N = 11, 11          # image dimensions
d1, d2 = 1, 1          # row and column delays
lam = 0.1              # EWMA smoothing parameter
L0 = 370.0             # target in-control ARL

ic = ICSTS(M, N, Normal(0, 1))   # in-control: spatial white noise
ICSTS(11, 11, Distributions.Normal{Float64}(μ=0.0, σ=1.0))

Computing the critical limit

cl = cl_sop(
    ic, lam, L0, 0.026, d1, d2;
    chart_choice=TauTilde(),
    reps_bracket=500, reps_final=2_000, bracket_step=0.002, seed=42,
    verbose=true
)

============================================================
 STEP 1: Bracket Search
 reps = 500 | truncation cap = 18500 | step size = 0.002
============================================================
  cl = 0.026 | ARL = 121.248
  cl = 0.028 | ARL = 179.01
  cl = 0.03 | ARL = 260.286
  cl = 0.032 | ARL = 414.918
------------------------------------------------------------
 Bracket found: [0.026, 0.032]
------------------------------------------------------------

============================================================
 STEP 2: Root Finding via ITP
 reps = 2000 | no truncation | tolerance = 0.0001
============================================================
  cl = 0.026 | ARL = 116.428
  cl = 0.032 | ARL = 393.329
  cl = 0.031487 | ARL = 354.9945
  cl = 0.031688 | ARL = 366.2585
  cl = 0.031731 | ARL = 368.107
  cl = 0.031751 | ARL = 370.267
Results of univariate zero finding:

* Converged to: 0.03175136751000527
* Algorithm: Roots.ITP{Float64, Int64}(0.2, 2, 1)
* iterations: 4
* function evaluations ≈ 6
* stopped as x_n ≈ x_{n-1} using atol=xatol, rtol=xrtol

Trace:
(a₀, b₀) = ( 0.025999999999999999, 0.032000000000000001 )
(a₁, b₁) = ( 0.031487298033593236, 0.032000000000000001 )
(a₂, b₂) = ( 0.031688040568479577, 0.032000000000000001 )
(a₃, b₃) = ( 0.031731176938555358, 0.032000000000000001 )
(a₄, b₄) = ( 0.031731176938555358, 0.03175136751000527 )

0.03175136751000527

verbose=true prints cl and the associated ARL at every evaluation, so both phases of the search can be followed: the bracketing steps away from cl_init until the ARL crosses L0, then the ITP algorithm refines the bracket. cl_init = 0.026 only has to be in the right ballpark, but it should sit below the root, and bracket_step should be of the order of the remaining distance — otherwise the bracketing crawls.

Increase reps_final for a more precise limit, and pass verbose=false (the default) to get the limit alone.

Verifying the in-control ARL

arl_sop_ic(ic, lam, cl, d1, d2, 5_000; chart_choice=TauTilde())
(376.8928, 5.247841269439258)

The returned tuple is (ARL, standard error); the ARL is within a few standard errors of L0.

Out-of-control performance

The same cl is used to monitor spatially dependent images. Note that the out-of-control DGP must generate images of the same size as the in-control process.

oc_dgps = [
    ("SAR(1,1), α = (0.1, 0.1, 0.0)", SAR11((0.1, 0.1, 0.0), M, N, Normal(0, 1), nothing, 20)),
    ("SAR(1,1), α = (0.2, 0.2, 0.0)", SAR11((0.2, 0.2, 0.0), M, N, Normal(0, 1), nothing, 20)),
    ("SAR(1,1), α = (0.3, 0.3, 0.1)", SAR11((0.3, 0.3, 0.1), M, N, Normal(0, 1), nothing, 20)),
    ("SQMA(1,1), β = (0.5, 0.3, 0.2)", SQMA11((0.5, 0.3, 0.2), (1, 1, 2), M, N, Normal(0, 1), nothing)),
]
oc_res = [
    arl_sop_oc(dgp, lam, cl, d1, d2, 5_000; chart_choice=TauTilde())
    for (_, dgp) in oc_dgps
]
DataFrame(
    Process=[n for (n, _) in oc_dgps],
    ARL=round.([r[1] for r in oc_res], digits=2),
    se=round.([r[2] for r in oc_res], digits=3),
)
4×3 DataFrame
Row Process ARL se
String Float64 Float64
1 SAR(1,1), α = (0.1, 0.1, 0.0) 17.51 0.142
2 SAR(1,1), α = (0.2, 0.2, 0.0) 6.98 0.036
3 SAR(1,1), α = (0.3, 0.3, 0.1) 5.29 0.022
4 SQMA(1,1), β = (0.5, 0.3, 0.2) 3.75 0.012

Even weak spatial dependence is detected within a few images, because every image contributes \((M - d_1)(N - d_2) = 100\) patterns.

Monitoring an observed image sequence

monitor_sop applied to a 3D array computes the EWMA chart statistics of an observed image sequence (the same ones stat_sop returns), applies the two sided alarm rule and records the first signal. The sequence below is spatial white noise up to \(t = 20\) and spatially dependent afterwards.

Random.seed!(42)
T, T_change = 40, 20
images = randn(M, N, T)

φ = 0.4
for t in (T_change+1):T, i in 2:M, j in 2:N
    images[i, j, t] = φ * (images[i-1, j, t] + images[i, j-1, t]) / 2 + randn()
end

res = monitor_sop(images, lam, cl, d1, d2; chart_choice=TauTilde())
ControlChartResult
  Family:           spatial ordinal patterns (SOP)
  Chart:            TauTilde()
  λ:                0.1
  Control limit:    0.0318
  Alarm rule:       |stat| > cl
  Statistics:       40 (images 1 to 40)
  First alarm:      image 23 (statistic -0.0392)

Each image contributes one statistic, so res.alarm is the index of the first image at which the chart signals:

(first_alarm=res.alarm, delay=res.alarm - T_change,
 stat_at_alarm=round(res.stats[res.alarm], digits=4))
(first_alarm = 23, delay = 3, stat_at_alarm = -0.0392)

The chart passes through the white-noise part without a false alarm and signals shortly after the change. With CairoMakie loaded, plot shows the statistics against the control limits, shades the rejection region and marks the first alarm. See the plotting tutorial for the options.

plot(res)

References

Adämmer, Philipp, Philipp Wittenberg, Christian H Weiß, and Murat Caner Testik. 2026. “Nonparametric Monitoring of Spatial Dependence.” Technometrics 68 (2): 267–81.