我假设你得到了一个长度为观察值的向量len正如您的示例中绘制的那样，并且您希望确定k跳跃和k跳跃尺寸。 （或者也许我误解了你；但你并没有真正说出你想要实现的目标。）
下面我将概述一个使用本地搜索的解决方案。我从您的示例数据开始：
x <- seq(0, 10, by = 0.01)
y <- staircase(x,
c(1,2,2,5),
c(2,5,2,1)) + rnorm(length(x), mean = 0, sd = 0.2)
解决方案是一个列表职位 and sizes的跳跃。请注意，我使用向量来存储这些数据，因为当你有 20 次跳跃时，定义变量会变得很麻烦。
示例（随机）解决方案：
k <- 5 ## number of jumps
len <- length(x)
sol <- list(position = sample(len, size = k),
size = runif(k))
## $position
## [1] 89 236 859 885 730
##
## $size
## [1] 0.2377453 0.2108495 0.3404345 0.4626004 0.6944078
我们需要一个目标函数来计算解决方案的质量。我还定义了一个简单的辅助函数stairs，由目标函数使用。
目标函数abs_diff计算拟合序列之间的平均绝对差（由解定义）和y.
stairs <- function(len, position, size) {
ans <- numeric(len)
ans[position] <- size
cumsum(ans)
}
abs_diff <- function(sol, y, stairs, ...) {
yy <- stairs(length(y), sol$position, sol$size)
sum(abs(y - yy))/length(y)
}
现在是本地搜索的关键组件：用于改进解决方案的邻域函数。邻域函数采用一个解决方案并对其进行轻微更改。在这里，它会选择一个position or a size并稍微修改一下。
neighbour <- function(sol, len, ...) {
p <- sol$position
s <- sol$size
if (runif(1) > 0.5) {
## either move one of the positions ...
i <- sample.int(length(p), size = 1)
p[i] <- p[i] + sample(-25:25, size = 1)
p[i] <- min(max(1, p[i]), len)
} else {
## ... or change a jump size
i <- sample.int(length(s), size = 1)
s[i] <- s[i] + runif(1, min = -s[i], max = 1)
}
list(position = p, size = s)
}
调用示例：此处新解决方案的第一个跳跃大小已更改。
## > sol
## $position
## [1] 89 236 859 885 730
##
## $size
## [1] 0.2377453 0.2108495 0.3404345 0.4626004 0.6944078
##
## > neighbour(sol, len)
## $position
## [1] 89 236 859 885 730
##
## $size
## [1] 0.2127044 0.2108495 0.3404345 0.4626004 0.6944078
我仍然负责本地搜索。
library("NMOF")
sol.ls <- LSopt(abs_diff,
list(x0 = sol, nI = 50000, neighbour = neighbour),
stairs = stairs,
len = len,
y = y)
我们可以绘制解决方案：拟合线显示为蓝色。
plot(x, y)
lines(x, stairs(len, sol.ls$xbest$position, sol.ls$xbest$size),
col = "blue", type = "S")
