在 R 中构建经验累积分布函数和数据插值答案

【问题标题】：Building an empirical cumulative distribution function and data interpolation in R在 R 中构建经验累积分布函数和数据插值
【发布时间】：2015-10-31 03:40:32
【问题描述】：

这是我正在使用的示例数据框

 level    Income    cumpop
 1      17995.50  0.028405
 2      20994.75  0.065550
 3      29992.50  0.876185
 4      41989.50  2.364170
 5      53986.50  4.267305
 6      65983.50  6.323390
 7      77980.51  8.357625
 8      89977.50 10.238910
 9     101974.50 11.923545
10     113971.51 13.389680
11     125968.49 14.659165
12     137965.50 15.753850
13     149962.52 16.673735
14     161959.50 17.438485
15     173956.50 18.093985
16     185953.52 18.640235
17     197950.52 19.099085
18     209947.52 19.514235
19     221944.50 19.863835
20     233941.50 20.169735
21     251936.98 20.628585
22     275931.00 20.936670
23     383904.00 21.850000

这个特定国家的全部人口已按收入分类，并分为 23 个相应的“级别”。 Income 变量是该级别所有成员的平均收入（这与说第 10 个百分位数的收入是 17995.50 很重要）。

但是每个级别的人口规模是不一致的（如果您查看cumpop 的差异，即累积人口，您会注意到这一点）。最终，我想构建一个 10 行数据框，为变量 Income 提供插值的十分位数，例如，我们可以说“最贫穷的 10% 人口平均有 28,000 人”或“人口中第 20% 到第 30% 的人平均赚 41,000 人”等等。所以我想有效地将这 23 个级别减少到 10 个人口规模相等的级别（以 cumpop[23] 作为总人口），这需要一些插值。

我环顾四周寻找一个库，它可以进行这种经验累积分布函数生成/插值，似乎ecdf 非常有用，但我不确定如何将它应用于Income 受@ 987654327@ 如上所述。

非常感谢这里的一些方向。

【问题讨论】：

标签： r interpolation ecdf

【解决方案1】：

使用黄土插值的快速而肮脏的解决方案。跨度设置得非常短，以确保基本上完美契合，遗憾的是，这也使任何错误术语变得毫无意义。尝试适当的回归可能值得。

incdist <- read.table("inc.txt", header=TRUE)

fit <- loess(incdist$Income~incdist$cumpop, span=0.2)
V2 <- predict(fit, seq(0, max(incdist$cumpop)*9/10, max(incdist$cumpop)/10))
V1 <- seq(0, max(incdist$cumpop)*9/10, max(incdist$cumpop)/10)
pred <- data.frame(V1, V2)

par(mar=c(5, 5.5, 4, 2) + 0.1)

plot(incdist$Income~incdist$cumpop, type="n", xaxt="n", yaxt="n",
    xlab="percentile", ylab=expression(frac("average income",1000)),
    main="income distribution")

abline(h=V2, v=V1[-1], col="grey")
points(incdist$Income~incdist$cumpop, col="grey")
lines(loess(incdist$Income~incdist$cumpop, span=0.2), col="red")
points(pred, col="blue", cex=1.5, pch=9)
axis(side=1, at=V1[-1], labels=c(1:9)*10)
axis(side=2, at=V2, labels=round(V2/1000), las=1)

【讨论】：