如何只拟合数据集的线性部分？

Question

p=(-50:50)^2
y=c(p, 2500+10*(1:99), p+1000)
plot(seq_along(y), y+100*rnorm(length(y)))

假设我有一个像上面这样的数据集，其中只有数据的一个子集是线性的。像lm()R 中的简单线性回归不能智能地找出适合线性拟合的区域（在本例中为 100 到 200）。

如何找出数据的哪一部分是线性的并仅在这个数据集子集中执行拟合？欢迎使用 R 和 python 中的解决方案。

请注意，上面显示的日期只是一个示例，只要它包含线性部分，该方法就应该对任意数据集具有鲁棒性。当有多个线性部分时，还应显示那些多个线性部分。如果没有线性部分，它应该显示没有找到线性部分。

编辑：一般来说，统计方法可能不适合稳健地解决这个问题。我添加了计算机视觉和机器学习标签。也许这些领域中的方法通常更适合稳健地解决这个问题？

Answer 1

我不知道有什么好的内置方法可以做到这一点，正如 Ben Bolker 和其他人指出的那样，这不是一个以稳健、可概括的方式回答的简单问题。也就是说，我使用蛮力方法在这个特定问题上取得了一些成功。因为我对tidyverse语法更熟悉，所以我使用了它，但我确信这可以在基础 R 中以类似的方式完成。

首先，我根据x序列的开始和长度创建了一个要探索的范围网格。根据您要执行的计算量调整粒度。对于快速方法，我使用每 5x和5length的倍数的 s。这给了我 1,830 个范围x，我附加了相关y的 。然后我将x和嵌套y到一个新列中data。

# From OP
p=(-50:50)^2
y=c(p, 2500+10*(1:99), p+1000)


library(tidyverse); library(broom)

df1 <- data.frame(x = seq_along(y), y = y+100*rnorm(length(y)))

df1_ranges = crossing(start  = seq.int(1, max(df1$x), by = 5), 
                      length = seq.int(5, 300, by = 5)) %>%
    mutate(end = start + length - 1) %>%
    filter(end <= max(df1$x)) %>%     # only keep ranges within the data
    uncount(length, .id = "x") %>%    # for each x, put in "length" many rows
    mutate(x = start + x - 1) %>%     # update x to run from "start" to "end"
    left_join(df1) %>%
    nest(data = c(x, y))

不是我可以lm对每个范围进行回归。这在我的电脑上大约需要 9 秒。您可以通过查看更少的不同范围或更聪明地了解搜索空间来加快速度。

df1_regressions <- df1_ranges %>%
    mutate(fit = map(data, ~lm(y~x, data = .x)),   # run lm's
           glance = map(fit, glance),              # summary of fit
           tidied = map(fit, tidy))                # extract coefficients

跳到后面，对于这个例子，具有最佳线性拟合的区域具有最低的回归项标准误差。果然，这确定了正确的位置，范围从大约 100 到 200。

df1_tidied <- df1_regressions %>%
    select(start:end, tidied) %>%
    unnest(tidied) %>%
    filter(term == "x")

df1_tidied %>%
    ggplot(aes(x =  start, y = end-start, fill = 1/std.error)) +
    geom_tile() +
    geom_text(data = . %>% filter(std.error == min(std.error)) %>% 
              mutate(text = glue::glue("({start}, {end-start})")), 
          aes(label = text), color = "white", vjust = -0.5) +
    scale_fill_viridis_c(direction = -1, option = "C")

哇！既然已经解决了，我们可以按照您最初的要求进行操作，并查看该部分的拟合回归。

df1_tidied %>% 
    slice_min(std.error) %>%
    select(start,end) %>%
    left_join(df1_ranges) %>%
    mutate(fit = map(data, ~lm(y~x, data = .x)),
           augment = map(fit, augment)) %>% 
    unnest(augment) -> df1_fitted

ggplot(df1, aes(x,y)) + 
    geom_point() +
    geom_line(data = df1_fitted, aes(y = .fitted), color = "red", size = 2)

Answer 2

dpseg在 dpseg 包中尝试一下。我们将最小长度限制为 50，以避免偶然发生的短线性拉伸。还有其他可用的调整参数。看?dpseg包装中附带的小插图。

为了使输入可重现，我们需要使用set.seed并在最后的注释中完成了这一点。

library(dpseg)
segs <- dpseg(x = x, y = y, minl = 50); segs
## ... this output is shown just before the image ...
subset(segs$segments, var < 20000)
##    x1  x2 start end intercept    slope        r2      var
## 3 116 203   116 203  1458.242 10.15865 0.8613225 10844.35

plot(segs)

给出以下结果，我们看到上面输出的第三段的方差最小。

> segs

calculating recursion for 301 datapoints

dynamic programming-based segmentation of 301 xy data points:

   x1  x2 start end  intercept     slope        r2      var
1   1  50     1  50   2165.902 -51.13574 0.9212552 47495.24
2  50 116    50 116  -2928.772  50.00892 0.9521128 47756.06
3 116 203   116 203   1458.242  10.15865 0.8613225 10844.35
4 203 252   203 252  12533.408 -47.39630 0.9189915 42079.16
5 252 301   252 301 -12405.806  51.67657 0.9261061 45278.70

Parameters: type: var; minl: 50; maxl: 301; P: 0; jumps: 0 

笔记

set.seed(123)
p <- (-50:50)^2
y <- c(p, 2500+10*(1:99), p+1000)
y <- y+100*rnorm(length(y))
x <- seq_along(y)