【问题标题】:How to only fit the linear portion of a dataset?如何只拟合数据集的线性部分?
【发布时间】:2021-10-27 19:42:33
【问题描述】:
p=(-50:50)^2
y=c(p, 2500+10*(1:99), p+1000)
plot(seq_along(y), y+100*rnorm(length(y)))

假设我有一个像上面这样的数据集,其中只有一部分数据是线性的。像 R 中的 lm() 这样的普通线性回归不能智能地找出适合线性拟合的区域(本例中为 100 到 200)。

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

注意,上面显示的日期只是一个示例,该方法对于任意数据集应该是稳健的,只要它包含线性部分。当有多个线性部分时,它也应该显示那些多个线性部分。如果没有线性部分,它应该显示没有找到线性部分。

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

【问题讨论】:

  • 看起来更像是一个统计问题,因为您明确询问的是一种方法而不是一种特定的编码方式
  • 示例中的线性区域在 100 到 200 之间。这部分的数据可以很好的与lm()吻合。
  • 不是一个容易以稳健/一般的方式回答的问题。
  • strucchange 包适合 linear 断点模型,但这已经是统计模型的广泛而深入的领域。推广到分段多项式模型(例如)将非常具有挑战性。
  • 我认为这只会对 OP 的满意度负责如果有人已经为此类问题编写了一个很好的通用求解器(我对此表示怀疑;正如我之前所说,一般来说,这是一个很难解决的问题)。我使用segmented 包尝试了一些拟合来拟合分段二次模型,这很快就表明它会很困难。

标签: machine-learning computer-vision linear-regression curve-fitting


【解决方案1】:

试试 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)

【讨论】:

  • 此方法无法识别 0->50, 50->100, 200->250, 250->300 是抛物线区域的一部分。它们在人眼中非常明显。所以这是朝着解决方案迈出的一大步,但它仍然不是一个理想的解决方案(我知道理想的解决方案可能很困难,这就是为什么在问题中添加计算机视觉标签的原因。)
  • 我认为很明显您会采用方差较小的分段,但显然并非如此,所以我明确添加了阈值 segs 的最后一步以给出线性分段 - 只是在这种情况下的第三段,根据需要。调整参数以适应。
  • 但是测试小变量是不够的。该点是抛物线区域中的线拟合,具有系统偏差,无法仅由 var 捕获。
  • 它适用于测试数据,如果您花一些时间在数据上设置调整参数,我认为您可以使用此框架来解决一般问题。非线性段可能具有更高的方差。请注意,dpseg 的 scoref 参数可用于提供替代评分,如果您不喜欢差异,可以花一些时间尝试不同的评分。
【解决方案2】:

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

首先,我根据起始x 和序列的长度创建了一个范围网格以供探索。根据您想要执行的计算量调整粒度。对于一个快速的方法,我使用了每 5 个 xlengths,它们是 5 的倍数。这给了我 1,830 个 x 范围,我在其中附加了相关的 y。然后我将xy 嵌套到一个新列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)

【讨论】:

  • 这是一个比应得的问题更好的答案,恕我直言。
  • 我认为这是一件有趣的事情。
  • 这实际上是一个很好的答案。
  • 我收到以下错误。如何解决? ``` R> df1_regressions mutate() 列 glance 有问题。 ℹglance = map(fit, glance)。 ✖ 找不到对象 'glance' 运行 rlang::last_error() 以查看发生错误的位置。 R> rlang::last_error() mutate()glance 有问题。 ℹglance = map(fit, glance)。 ✖ 找不到对象'glance' 回溯:1.%&gt;%(...) 9.base::.handleSimpleError(...) 10.dplyr:::h(simpleError(msg, call)) 运行rlang::last_trace()查看完整的上下文。 ```
  • tidyverse 有很多缺点。 stackoverflow.com/questions/61393345/…所以我不熟悉,不使用。特别是,它太冗长了。等价基数R应该更简洁吧?欢迎其他可以提供等效基础 R 解决方案的人。
猜你喜欢
  • 2020-10-12
  • 1970-01-01
  • 2016-03-09
  • 2019-03-15
  • 2010-10-09
  • 1970-01-01
  • 2018-12-18
  • 2020-07-27
  • 1970-01-01
相关资源
最近更新 更多