将指数分布拟合到频率 table

Question

我有以下数据集：

intervals <- c("0-10", "10-20", "20-30", "30-40", "40-50", "50-75", "75-100", ">100")
int.mean <- c(5.5, 14.3, 24.9, 35.4, 45.2, 63.1, 86.1, NA)
freq <- c(165, 90, 55, 25, 20, 35, 30, 15)

data <- data.frame(intervals, int.mean, freq)

我希望对数据进行指数分布拟合，以预测值超过 150 的概率并具有一定的置信度。我可以按如下方式拟合分布：

library(MASS)
fittedexp <- fitdistr(na.exclude(data$int.mean), "exponential")

然而，这并没有考虑到频率，所以我不确定我这样做是否正确。然后我计划使用 optim 函数来创建估计概率的置信区间。

Answer 1

您可以使用 freq 变量扩展数据，然后拟合分布

data.expand <- data[rep(seq_len(nrow(data)), times=data$freq), ]
head(data.expand, 3); tail(data.expand, 3)

    intervals int.mean freq             intervals int.mean freq
1        0-10      5.5  165        8.12      >100       NA   15
1.1      0-10      5.5  165        8.13      >100       NA   15
1.2      0-10      5.5  165        8.14      >100       NA   15

library(MASS)
with(subset(data.expand, subset=!is.na(int.mean)),
        fitdistr(int.mean,densfun="exponential")
)

      rate    
  0.041401745 
 (0.002020198)

Answer 2

您正在处理一个分类变量 "intervals"，它根据您从中获取断点的假定基础连续变量创建离散的计数观察。有点混乱的数据情况。从技术上讲，您有 interval-censored 数据。 ~~但是，如果您将指数分布作为假设，那么您计算的那些 "means" 实际上是中点，但它们不会是指数分布变量的均值。~~ 见下文对于我对 int.means 观察结果的修改意见。（所以现在我将扩展我的原始评论以包含一些 R 代码。）

如果我们将间隔的端点作为中断变量，并计算观察数据中的比例，我们有：

 brks <- c(0, 10,20,30,40,50,75,100,Inf)
 freq <- c(165, 90, 55, 25, 20, 35, 30, 15)
 prop<- freq/sum(freq)
 prop
#-----
[1] 0.37931034 0.20689655 0.12643678 0.05747126 0.04597701 0.08045977 0.06896552 0.03448276
round(prop,2)
[1] 0.38 0.21 0.13 0.06 0.05 0.08 0.07 0.03

然后我们可以显示具有相似均值的指数分布变量可能"look like"（就比例而言）如果分到这些区间：

 table( findInterval( rexp(100, 1/15), brks) )/100

   1    2    3    4    5    6    7 
0.47 0.24 0.12 0.08 0.04 0.04 0.01

所以我们可能想尝试一个高于 15 的平均值，比如说 20？

> table( findInterval( rexp(100, 1/20), brks) )/100

   1    2    3    4    5    6    7    8 
0.37 0.24 0.13 0.09 0.07 0.07 0.02 0.01 
> round(prop,2)
[1] 0.38 0.21 0.13 0.06 0.05 0.08 0.07 0.03

所以您可以很好地拟合观测值的低端，但指数分布的变量似乎有点 "thinner" 尾巴。由于您对数据的高端感兴趣，因此您可能希望在高端获得更好的拟合，但这会扰乱您的统计原则置信区间目标。你有点卡住了，因为你的数据不是一组正确的 "exponential" 观察值。（将模拟大小增加到 1000 以减少噪声的影响。）

> table( findInterval( rexp(1000, 1/25), brks) )/1000

    1     2     3     4     5     6     7     8 
0.329 0.222 0.141 0.103 0.056 0.094 0.034 0.021 
> round(prop,2)
[1] 0.38 0.21 0.13 0.06 0.05 0.08 0.07 0.03

那里的身材看起来并不糟糕。如果指数分布的速率参数是 1/25，那么这将是大于 150 的观测值的比例：

 1-pexp(150, 1/25)
#[1] 0.002478752

可能有用：http://jsdajournal.springeropen.com/articles/10.1186/s40488-015-0028-6

您也可以尝试在 CrossValidated.com 上搜索存在一些先前讨论的地方。

编辑：我最初认为那些 int.means 值是间隔边界的中点，但显然不是这样，因为它们似乎接近中点但有很大的抖动围绕中点。此外，这些值与指数分布不一致，因为在人口最多的区间 (0-10) 中，观测值应该在中点的 "left" 处，甚至不在中点的左侧。它可能应该是 4.0 或 4.5，但肯定不会高到 5.5。它表明一些其他分布是这个物理过程的基础，也许是某种 Gamma 分布，它会在零附近下降到零，但在 0-10 间隔的早期达到峰值，然后有一个更长的尾巴。

将指数分布拟合到频率 table

Fitting exponential distribution to frequency table

r

exponential-distribution