8

我正在搜索 John Tukey 算法，该算法在我使用 R 的线性回归上计算“阻力线”或“中位数线”。

邮件列表上的一个学生用这些术语解释了这个算法：

“它的计算方式是将数据分成三组，找到每组的x-median和y-median值（称为汇总点），然后用这三个汇总点确定线。最外面的两个汇总点决定斜率，所有点的平均值决定截距。”

关于 John tukey 的好奇中位数的文章：http: //www.johndcook.com/blog/2009/06/23/tukey-median-ninther/

你知道我在哪里可以找到这个算法或 R 函数吗？在哪个包中，非常感谢！

4

3 回答 3

12

这里有关于如何计算中位数线的描述。一个R实现是

median_median_line <- function(x, y, data)
{
  if(!missing(data))
  {
    x <- eval(substitute(x), data) 
    y <- eval(substitute(y), data) 
  }
  
  stopifnot(length(x) == length(y))

  #Step 1
  one_third_length <- floor(length(x) / 3)
  groups <- rep(1:3, times = switch((length(x) %% 3) + 1,
     one_third_length,
     c(one_third_length, one_third_length + 1, one_third_length),
     c(one_third_length + 1, one_third_length, one_third_length + 1)
  ))

  #Step 2
  x <- sort(x)
  y <- sort(y)
  
  #Step 3
  median_x <- tapply(x, groups, median)                                 
  median_y <- tapply(y, groups, median)

  #Step 4
  slope <- (median_y[3] - median_y[1]) / (median_x[3] - median_x[1])
  intercept <- median_y[1] - slope * median_x[1]

  #Step 5
  middle_prediction <- intercept + slope * median_x[2]
  intercept <- intercept + (median_y[2] - middle_prediction) / 3
  c(intercept = unname(intercept), slope = unname(slope))
}

为了测试它，这里有一个例子：

dfr <- data.frame(
  time = c(.16, .24, .25, .30, .30, .32, .36, .36, .50, .50, .57, .61, .61, .68, .72, .72, .83, .88, .89),
  distance = c(12.1, 29.8, 32.7, 42.8, 44.2, 55.8, 63.5, 65.1, 124.6, 129.7, 150.2, 182.2, 189.4, 220.4, 250.4, 261.0, 334.5, 375.5, 399.1))
  
median_median_line(time, distance, dfr) 
#intercept     slope 
#   -113.6     520.0

请注意指定组的稍微奇怪的方式。这些说明对如何定义组大小非常挑剔，因此更明显的方法cut(x, quantile(x, seq.int(0, 1, 1/3)))不起作用。

于 2010-07-12T13:36:28.903 回答
2

我参加聚会有点晚了，但是您尝试过 stats 包中的 line() 吗？

从帮助文件：

价值

“tukeyline”类的对象。

参考

图基，JW（1977）。探索性数据分析，阅读马萨诸塞州：Addison-Wesley。

于 2013-05-31T09:16:02.043 回答
2

作为 R Core 团队的一员，我现在已经深入研究了源代码，也研究了它的历史。

结论：在 19961997 年添加的源 C 源代码，当时 R 仍被称为 alpha（和 0.14alpha 版本左右）已经计算了不太正确的分位数......对于某些样本大小。

更多关于 R 邮件列表的信息（还没有）。

于 2014-12-19T17:18:37.587 回答