r - 识别数据中的几个精确线性关系

Question

我有一些二维数据，其中大量的行正好服从少数线性关系中的一种。绘制数据时很容易识别线条:

如何识别这些线的斜率和截距？

虽然应用哪种线性关系是一个基于另一个变量的确定性过程，但那个变量已经丢失了。我不在乎我无法预测新值；我只想要所有的斜率和截距。

如果截距为零，则算法相对容易。只需为每个点计算 r = y/x，将其四舍五入到一定精度，然后确定最频繁的 r。但是，当截距不为零时，这不会推广。

Reproducible data:

library(data.table)
div <- function(i, d) {
  {i %% d} == 0L
}

DT <- data.table(x = runif(1e6, 1, 100e3), i = seq_len(1e6))
DT[, y := 0.8 * x + 23333]
DT[div(i, 3), y := 0.3 * x + 14444]
DT[div(i, 7), y := 1.7 * x + 8888]
DT[1:50e3, y := y + runif(.N, -20e3, 20e3)]

我尝试做的一个过程是执行交叉连接，计算点样本与所有其他点之间的斜率。在这种情况下，它确实识别了斜坡；然而，这只需要少数点脱线并且可能有点低效。

CJ1 <- function(seq., siz = 500) {
  CJ(i1 = seq.,
     i2 = sample.int(1e6, size = siz)) %>%
    .[DT, on = "i1==i", nomatch = 0L] %>%
    .[DT[, .(x1 = x, y1 = y, i2 = i)], on = "i2", nomatch = 0L] %>%
    .[, m := round((y - y1) / (x - x1), 3)] %>%
    .[, .N, keyby = .(m)] %>%
    .[order(-N)] %>%
    .[N > (2 * N[20])]
}

是否有任何成熟的建模技术来提取这种线性关系？

Answer 1

我将示例数据框的行数减少到 1000，如下所示，因为这与您提供的绘图示例更相似。

library(data.table)

div <- function(i, d) {
  {i %% d} == 0L
}

set.seed(1)

DT <- data.table(x = runif(1e3, 1, 100e3), i = seq_len(1e3))
DT[, y := 0.8 * x + 23333]
DT[div(i, 3), y := 0.3 * x + 14444]
DT[div(i, 7), y := 1.7 * x + 8888]
DT[1:(1e2), y := y + runif(.N, -20e3, 20e3)]

library(ggplot2)

# See the scatter plot of x and y
ggplot(DT, aes(x = x, y = y)) +
  geom_point()

如您所述，如果直线穿过原点。斜率可以计算为 y/x。因此，我们可以遍历一系列可能的截距值，例如 1:30000，计算 y 减去潜在截距，然后计算所有行的斜率。如果真实截距显示，数据集中将有一个斜率变得丰富。我们可以使用一个函数来找到模式，从而我们可以确定截距和斜率。下面是一个代码示例。

# A function to find the mode
# From https://stackoverflow.com/a/8189441/7669809
Modes <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux))
  ux[tab == max(tab)]
}

# A for loop to loop through a series of intercept and determine the number of mode
answer <- numeric()

for (i in 1:30000){
  y_temp <- DT$y - i
  r <- y_temp/DT$x
  M <- Modes(r)
  if (length(M) == 1){
    answer[[as.character(i)]] <- M
  }
}

最终输出是一个以斜率作为元素并以截距作为关联名称的向量。

answer
# 8888 14444 23333 
#  1.7   0.3   0.8