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)))
# build adjacency graph
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)
