Choosing mean vs median for noisy sensor rollups
A tri-axial accelerometer bolted to a gearbox on a stamping press reports RMS vibration velocity every 200 ms, and most of the time that signal is quiet and boring. Then the tool strikes: a 40 ms impact transient spikes the reading twenty times above baseline, the press indexes, and the channel settles again. When a 1-minute rollup collapses those ~300 raw points with fn: mean, a single impact drags the aggregate off-scale, your condition-monitoring dashboard lights up amber, and the on-call engineer walks the floor to find a machine running perfectly. The problem is not the threshold — it is the aggregation function feeding the threshold. Choosing between mean, median, and a trimmed quantile is a decision about which statistic survives contact with outliers, and it is the first knob to turn before you touch alert levels in Threshold Tuning for Aggregation. This page shows how to quantify the skew, pick the robust statistic in Flux, and keep multi-tier rollups honest.
Prerequisites
Solution walkthrough
Step 1 — Quantify how far the mean is being pulled
Do not switch functions on a hunch. Measure the divergence between mean and median per window so the choice is evidence-driven. The gap between the two is a direct read of skew: a symmetric, quiet window has mean ≈ median, while a window containing an impact shows the mean sitting well above the median.
import "join"
raw = from(bucket: "vibration-raw")
|> range(start: -6h)
|> filter(fn: (r) => r._measurement == "vibration" and r._field == "velocity_mms")
|> filter(fn: (r) => r.machine_id == "press-07")
means = raw |> aggregateWindow(every: 1m, fn: mean, createEmpty: false)
meds = raw |> aggregateWindow(every: 1m, fn: median, createEmpty: false)
join.left(
left: means |> keep(columns: ["_time", "_value"]),
right: meds |> keep(columns: ["_time", "_value"]),
on: (l, r) => l._time == r._time,
as: (l, r) => ({_time: l._time, mean: l._value, median: r._value, skew: l._value - r._value}),
)
|> sort(columns: ["skew"], desc: true)
|> limit(n: 20)
Sorting by skew descending surfaces exactly the windows where the mean is being levered away from the bulk of the data. If the top rows show a skew of several mm/s against a baseline near 2 mm/s, a mean-based threshold is firing on transients, not sustained faults. Formally, one spike of magnitude $x_k$ inside an $n$-point window moves the mean by
$$ \Delta \bar{x} = \frac{x_k - \bar{x}_{\setminus k}}{n} $$
so with $n \approx 300$ raw points per minute a single 60 mm/s impact against a 2 mm/s floor still lifts the reported mean by roughly 0.19 mm/s — small per point, but a train of impacts compounds, and shorter windows amplify it. The median, by contrast, does not move at all until outliers occupy half the window.
Step 2 — Pick the robust statistic in the rollup task
Once you have decided the spikes are contamination for the trend channel, roll up with median (the 50th percentile) instead of mean. The task below writes a robust central-tendency series into the aggregate bucket that your threshold check then reads.
option task = {name: "rollup-vibration-1m", every: 1m, offset: 15s}
from(bucket: "vibration-raw")
|> range(start: -task.every)
|> filter(fn: (r) => r._measurement == "vibration" and r._field == "velocity_mms")
|> aggregateWindow(every: 1m, fn: median, createEmpty: false)
|> set(key: "_field", value: "velocity_p50")
|> to(bucket: "vibration-1m", org: "reliability-eng")
fn: median is shorthand for quantile(q: 0.5) and uses the approximate t-digest estimator by default — cheap and more than accurate enough for a trend line. If you want to keep the peak as a separate safety signal rather than discard it, run a second pass with fn: max into a velocity_peak field; that way the median tracks machine health while a genuine over-limit impact still trips its own hard threshold. The offset: 15s lets the last stragglers of each minute land before the window closes, the same late-data discipline that governs the wider downsampling aggregation pipeline. For contamination you want gone entirely — dead-sensor zeros, wiring glitches — reach instead for a trimmed reading:
|> aggregateWindow(
every: 1m,
fn: (column, tables=<-) => tables |> quantile(q: 0.9, column: column, method: "exact_selector"),
createEmpty: false,
)
A high quantile like the p90 rejects the bottom of the distribution and the single top-most spike alike, giving a “typical loaded” value. method: "exact_selector" returns an actual data point (useful when you need a real measured reading), whereas the default estimate_tdigest interpolates and is markedly faster on large windows — a trade-off that matters at fleet scale.
Step 3 — Keep multi-tier rollups additive
Here is the trap that bites weeks later: the median is not composable. A daily rollup cannot be built from 1,440 one-minute medians, because the median of medians is not the median of the raw data. The mean has the same limitation unless you carry weights. So when you cascade rollups, decide per field:
- Additive fields (sum, count): roll the daily tier from the hourly tier freely — sums compose.
- Mean central tendency: carry both
sumandcountup each tier and divide only at read time, so the day’s mean is $\left(\sum \text{sum}\right) / \left(\sum \text{count}\right)$, a true count-weighted mean rather than a mean-of-means. - Median / quantiles: recompute each tier directly from the raw (or a t-digest sketch), never from the tier below. In practice this means the daily median task reads
vibration-raw, notvibration-1m.
// daily robust rollup reads RAW, not the 1-minute tier
option task = {name: "rollup-vibration-1d", every: 1d, offset: 20m}
from(bucket: "vibration-raw")
|> range(start: -task.every)
|> filter(fn: (r) => r._measurement == "vibration" and r._field == "velocity_mms")
|> aggregateWindow(every: 1d, fn: median, createEmpty: false)
|> set(key: "_field", value: "velocity_p50")
|> to(bucket: "vibration-1d", org: "reliability-eng")
Reading raw for every tier is more expensive, which is why the choice of statistic and the choice of schedule are linked; the ordering and offset discipline that keeps cascaded tiers correct is covered in multi-tier rollup scheduling. If raw is already expired by the time the daily task runs, you must either widen raw retention or accept an approximate quantile stored as a t-digest.
Gotchas and edge cases
Median-of-medians silently understates variability. The most common mistake is pointing the daily task at the 1-minute rollup because it is cheaper. Each minute’s median has already smoothed away its own spikes, so the day’s median-of-medians is even flatter than reality — you lose the very transients a reliability team needs to see, with no error to warn you. Recompute robust statistics from raw, or store a mergeable t-digest sketch per window instead of a scalar.
median on t-digest can disagree with an exact selector at boundaries. The default estimate_tdigest method interpolates and may return a value that never appears in the data. That is fine for a trend line but wrong if a downstream job assumes the rollup is a real measured sample (for example, matching it back to a raw timestamp). Use method: "exact_selector" when you need a genuine data point, and budget for the extra memory it holds per window.
Discarding spikes can hide the fault you were hired to catch. Robustness cuts both ways. If those impacts are the failure mode — a cracked bearing throwing periodic knocks — a median rollup will happily report a healthy machine right up to seizure. When spikes are signal, keep a parallel max or p99 channel with its own hard threshold so the smoothed trend and the peak detector run side by side; never let one function serve both purposes.
Verification
Confirm the robust rollup tracks the bulk of the signal while the mean channel chases the spikes. Write both, then compare a known impact window:
union(tables: [
from(bucket: "vibration-1m")
|> range(start: -1h)
|> filter(fn: (r) => r._field == "velocity_p50"),
from(bucket: "vibration-raw")
|> range(start: -1h)
|> filter(fn: (r) => r._field == "velocity_mms")
|> aggregateWindow(every: 1m, fn: mean, createEmpty: false)
|> set(key: "_field", value: "velocity_mean"),
])
|> pivot(rowKey: ["_time"], columnKey: ["_field"], valueColumn: "_value")
|> map(fn: (r) => ({r with divergence: r.velocity_mean - r.velocity_p50}))
|> filter(fn: (r) => r.divergence > 5.0)
|> yield(name: "windows_where_mean_lies")
Any rows returned are minutes where the mean-based threshold would have fired but the median stayed calm — precisely the false alarms you set out to eliminate. A healthy result is a short list confined to real impact events, confirming the median rollup is reporting sustained machine condition rather than transients.
Related
- Threshold Tuning for Aggregation — how alert levels and hysteresis sit on top of the aggregate statistic chosen here.
- Multi-Tier Rollup Scheduling — cascading hourly and daily tiers without recomputing non-additive statistics incorrectly.
- Downsampling & Aggregation Pipeline Design — the end-to-end rollup architecture these tasks plug into.
Up one level: Threshold Tuning for Aggregation