r - 带有 2 个断点的断棒（或分段）回归

Question

我想用下一个数据估计一个函数的两个断点：

    df = data.frame (x = 1:180,
                y = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 2, 2, 4, 2, 2, 3, 2, 1, 2,0, 1, 0, 1, 4, 0, 1, 2, 3, 1, 1, 1, 0, 2, 0, 3,  2, 1, 1, 1, 1, 5, 4, 2, 1, 0, 2, 1, 1, 2, 0, 0, 2, 2, 1, 1, 1, 0, 0, 0, 0, 
                    2, 3, 0, 3, 2, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
# plotting y ~ x 
plot(df)

我知道该函数有两个断点，例如：

y = y1 if x < b1;
y = y2 if b1 < x < b2;
y = y3 if b2 < x;

我想找到b1并b2拟合一种具有以下形式的矩形函数

任何人都可以帮助我或指出我正确的方向吗？谢谢！

标签： rlinear-regressionpiecewise

Answer 1

1）kmeans试试kmeans这样：

set.seed(123)
km <- kmeans(df, 3, nstart = 25)

> fitted(km, "classes") # or equivalently km$cluster
  [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [38] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[149] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

> unique(fitted(km, "centers")) # or except for order km$centers
      x         y
3  30.5 0.5166667
1  90.5 0.9000000
2 150.5 0.0000000

> # groups are x = 1-60, 61-120 and 121-180
> simplify2array(tapply(df$x, km$cluster, range))
       1   2  3
[1,]  61 121  1
[2,] 120 180 60

plot(df, col = km$cluster)
lines(fitted(km)[, "y"] ~ x, df)

2）蛮力 另一种方法是蛮力方法，我们计算每对可能的断点，并选择线性模型中平方和最小的对。

grid <- subset(expand.grid(b1 = 1:180, b2 = 1:80), b1 < b2)

# the groups are [1, b1], (b1, b2], (b2, Inf)
fit <- function(b1, b2, x, y) {
   grp <- factor((x > b1) + (x > b2))
   lm(y ~ grp)
}

dv <- function(...) deviance(fit(...))

wx <- which.min(mapply(dv, grid$b1, grid$b2, MoreArgs = df))

grid[wx, ]
##       b1 b2
## 14264 44 80

plot(df)
lines(fitted(fit(grid$b1[wx], grid$b2[wx], x, y)) ~ x, df)