我有一个有点嘈杂的时间序列，它徘徊在不同的水平上。
例如以下数据：
我有可用的实线数据，我想获得虚线的估计值。它应该是分段常数。
什么算法适合在这里尝试？
到目前为止，我的想法徘徊在 0 度 P 样条曲线（但如何找出在哪里放置结？）或结构断裂模型。回归树是我目前最好的想法，但理想情况下，我会寻找一种方法，考虑到 y=250 处的两个水平具有相同的 y 值这一事实。如果我理解正确，回归树会将这两个区间分成两个不同的组，每个组都有不同的平均值。
生成它的 R 代码是这样的：
set.seed(20181118)
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
x = 1:400
y = true_fct(x) + rt(length(x), df=1)
plot(x, y, type="l")
lines(x, true_fct(x), lty=2, lwd=3)
处理此类噪声的一种简单、稳健的方法是计算中值。
短窗口上的滚动中位数将检测除最小跳跃之外的所有跳跃，而检测到的跳跃之间的间隔内的响应中位数将有力地估计其水平。（您可以用不受异常值影响的任何稳健估计替换后一个估计。）
您应该使用真实或模拟数据调整此方法，以达到可接受的错误率。例如，对于问题中的模拟，我发现使用第二个和第 98 个百分位数来设置检测跳跃的阈值是很好的。在其他情况下——例如当可能发生许多跳跃时——更多的中心百分位数会更好。
以下是显示 (a) 三个跳跃为红点和 (b) 四个估计水平为浅蓝色线的结果。
跳跃估计发生在索引 100、200、250（这正是模拟使它们发生的位置），结果水平估计在 199.6、249.8、300.0 和 250.2：均在真实基础值的 0.4 范围内。
这种出色的行为在重复模拟中仍然存在（set.seed在开始时删除命令）。
这是R代码。
#
# Rolling medians.
#
rollmed <- function(x, k=3) {
n <- length(x)
x.med <- sapply(1:(n-k+10), function(i) median(x[i + 0:(k-1)]))
l <- floor(k/2)
c(rep(NA, l), x.med, rep(NA, k-l))
}
y.med <- rollmed(y, k=5)
#
# Changepoint analysis.
#
dy <- diff(y.med)
fourths <- quantile(dy, c(1,49)/50, na.rm=TRUE)
thresholds <- fourths + diff(fourths)*2.5*c(-1,1)
jumps <- which(dy < thresholds[1] | dy > thresholds[2]) + 1
points(jumps, y.med[jumps], pch=21, bg="Red")
#
# Plotting.
#
limits <- c(1, jumps, length(y)+1)
y.hat <- rep(NA, length(jumps)+1)
for (i in 1:(length(jumps)+1)) {
j0 <- limits[i]
j1 <- limits[i+1]-1
y.hat[i] <- median(y[j0:j1])
lines(x[j0:j1], rep(y.hat[i], j1-j0+1), col="skyblue", lwd=2)
}
如果您仍然对使用 L0 惩罚进行平滑处理感兴趣，我会查看以下参考资料：“使用 L0 惩罚进行分段平滑的基因组变化可视化”-DOI：10.1371/journal.pone.0038230（对Whittaker smoother 可以在 P. Eilers 论文“A perfect smoother”中找到 - DOI: 10.1021/ac034173t）。当然，为了实现您的目标，您必须围绕该方法进行一些工作。
原则上，您需要 3 种成分：
当然，您还需要一种选择最佳平滑量的方法。这个例子是由我的木匠眼睛完成的。您可以使用 DOI 中的标准：10.1371/journal.pone.0038230（第 5 页，但我没有在您的示例中尝试过）。
您将在下面找到一个小代码。我留下了一些评论作为指导。
# Cross Validated example
rm(list = ls()); graphics.off(); cat("\014")
library(splines)
library(Matrix)
library(quantreg)
# The data
set.seed(20181118)
n = 400
x = 1:n
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
y = true_fct(x) + rt(length(x), df = 1)
# Prepare bases - Identity matrix (Whittaker)
# Can be changed for B-splines
B = diag(1, n, n)
# Prepare penalty - lambda parameter fix
nb = ncol(B)
D = diff(diag(1, nb, nb), diff = 1)
lambda = 1e2
# Solve standard Whittaker - for initial values
a = solve(t(B) %*% B + crossprod(D), t(B) %*% y, tol = 1e-50)
# est. loop with L0-Diff penalty as in DOI: 10.1371/journal.pone.0038230
p = 1e-6
nit = 100
beta = 1e-5
for (it in 1:nit) {
ao = a
# Penalty weights
w = (c(D %*% a) ^ 2 + beta ^ 2) ^ ((p - 2)/2)
W = diag(c(w))
# Matrix augmentation
cD = lambda * sqrt(W) %*% D
Bp = rbind(B, cD)
yp = c(y, 1:nrow(cD)*0)
# Update coefficients - rq.fit from quantreg
a = rq.fit(Bp, yp, tau = 0.5)$coef
# Check convergence and update
da = max(abs((a - ao)/ao))
cat(it, da, '\n')
if (da < 1e-6) break
}
# Fit
v = B %*% a
# Show results
plot(x, y, pch = 16, cex = 0.5)
lines(x, y, col = 8, lwd = 0.5)
lines(x, v, col = 'blue', lwd = 2)
lines(x, true_fct(x), col = 'red', lty = 2, lwd = 2)
legend("topright", legend = c("True Signal", "Smoothed signal"),
col = c("red", "blue"), lty = c(2, 1))
PS。这是我对交叉验证的第一个答案。我希望它足够有用和清晰:-)
我会考虑使用Ruey Tsay 的论文 Outliers, level shifts, and variance changes in time series Differencing model with AR1 and 21 outliers。
我们关闭了differencng，并专门调用了电平转换。
