我不知道有什么好的内置方法可以做到这一点,正如 Ben Bolker 和其他人所指出的那样,这不是一个以稳健、可概括的方式回答的直截了当的问题。也就是说,我使用蛮力方法在这个特定问题上取得了一些成功。因为我更习惯tidyverse 语法,所以我使用了它,但我确信这可以在基本 R 中以类似的方式完成。
首先,我根据起始x 和序列的长度创建了一个范围网格以供探索。根据您想要执行的计算量调整粒度。对于一个快速的方法,我使用了每 5 个 x 和 lengths,它们是 5 的倍数。这给了我 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)