# Detect constant (zero-slope) sections in a noisy step function

I have data that has the shape of a step function, which is obstructed by noise. I would like to identify the sections with constant slope in it.

The noise does not necessarily have 0 mean, and the steps are not always abrupt. Also, the points on the x-axis are not equidistantly spaced.

I was thinking about trying one-sided differentiators (like here) but not sure if they would work with non-equidistantly spaced data.

I feel like some sort of a global approach based on optimization would be an elegant way to solve that but so far I couldn't come up with one.

• Is it literally a step function, like it only produces certain discrete values? – David Marx Jan 18 '18 at 10:54
• I've added example plots. Wanted to add them initially but forgot. – Konstantin Jan 18 '18 at 13:01
• In the second plot, around 1900000, there is a really steep small step and a big, sloped step, would you count this as one or two steps? – Jan van der Vegt Jan 18 '18 at 13:57
• That's exactly one of the challenges, and also the reason why the differentiator hasn't worked very well. Those shall count as one step, as after the first steep small step the new level is not kept. Basically, this effect is due to the slow reaction of the sensor and the delay in the sensor value updates. – Konstantin Jan 18 '18 at 14:35

For what you're trying to do, first-order differences without interpolation should work just fine. Once you've done that, the problem reduces to a simple anomaly detection task. Here's a demonstration using topoligical anomaly detection. Briefly:

The technique is essentially a density based outlier detection algorithm that, instead of calculating local densities, constructs a graph of the data using nearest-neighbors. The algorithm is different from other kNN outlier detection algorithms in that instead of setting ‘k’ as a parameter, you instead set a maximal inter-observation distance (called the graph “resolution” by Gartley and Basener). If the distance between two points is less than the graph resolution, add an edge between those two observations to the graph. Once the full graph is constructed, determine which connected components comprise the “background” of the data by setting some threshold percentage of observations ‘p’: any components with fewer than ‘p’ observations is considered an anomalous component, and all the observations (nodes) in this component are outliers.

First, let's generate some data:

n_obs = 100
spike_perc = .1

###################

set.seed(123)

n_spikes = floor(spike_perc * n_obs)
inter_arrival = rexp(n_obs)
xv = cumsum(inter_arrival)

mu_spike = 5
sd_spike = sqrt(mu_spike)
sd_noise = mu_spike/75
spike_events = sample(n_obs, n_spikes)

flat_noise = rnorm(n_obs, 0, sd_noise)
spike_ampl = abs(rnorm(n_spikes, mu_spike, sd_spike))

yv = rep(0, n_obs)
yv[spike_events] = yv[spike_events]  + spike_ampl
yv = cumsum(yv)
yv = yv + flat_noise


Here's what our data looks like (spike events in red):

Calculate first-order differences (appending 0 to preserve index mapping):

dat = c(0,diff(yv))


Find us some outliers:

library(igraph)

# Helper function. Equivalently, use mefa::stack
row_col_from_condensed_index = function(n, ix){
nr=ceiling(n-(1+sqrt(1+4*(n^2-n-2*ix)))/2)
nc=n-(2*n-nr+1)*nr/2+ix+nr
cbind(nr,nc)
}

# Algorithm parameters, may require tuning
rq = .2 # distance quantile for thresholding (higher->fewer outliers)
p = .1  # percent of population necessary for a component to be flagged as an inlier
method='euclidean'

# calculate inter-observation distances and threshold to define adjacency matrix
n=NROW(dat)
d = dist(dat, method=method)
r = quantile(d, rq)
edges = t(sapply(which(d<=r), function(x) row_col_from_condensed_index(n,x)))

g = graph.empty(n, directed=FALSE)
g[from=edges[,1], to=edges[,2]]=1
V(g)$name = 1:n # Flag outliers based on component size components = clusters(g, mode="strong") csize = components$csize[components$membership] outliers = V(g)[csize<=p*n]$name


And here's how we did (flagged outliers in red):

If you go this route, you may need to tweak the rq and p parameters a bit. If you'd rather not deal with the quantile thresholding parameter (rq), an alternative approach is to use a nearest neighbor graph, in which case you'll have to tinker with the 'k' parameter (of KNN). That approach would replace the distance calculations and graph instantiation above with:

library(cccd)

g = nng(data.frame(dat), k=5, use.fnn=TRUE)

• I was stuck with the idea to apply differentiation, which gave me headache because if the noise. Obviously, using differences is the way to go here. Special thanks for pointing out TAD. – Konstantin Jan 19 '18 at 9:33