将指数分布拟合到频率表答案

【问题标题】：Fitting exponential distribution to frequency table将指数分布拟合到频率表
【发布时间】：2020-02-27 00:32:47
【问题描述】：

我有以下数据集：

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 函数来创建估计概率的置信区间。

【问题讨论】：

标签： r exponential-distribution

【解决方案1】：

您正在处理一个分类变量“区间”，它根据假定的基础连续变量创建一个离散的计数观察，您从中获取断点。那种凌乱的数据情况。从技术上讲，您有区间删失数据。 ~~但是，如果您有指数分布作为假设，那么您计算的那些“均值”实际上是中点，但它们不会被期望是指数分布变量的均值。~~请参阅下面的修改后的 cmets 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

然后，我们可以展示一个具有相似均值的指数分布变量如果合并到这些区间中可能“看起来”（就比例而言）：

 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

因此，您可以很好地拟合观测值的低端，但指数分布的变量似乎有一条“更细”的尾巴。由于您对数据的高端感兴趣，因此您可能希望在高端获得更好的拟合，但这会与您的统计原则置信区间目标相混淆。你有点卡住了，因为你的数据不是一组正确的“指数”观察。（将模拟大小增加到 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) 中，观察结果应该在中点的“左侧”，甚至不在中点的左侧。它可能应该是 4.0 或 4.5，但肯定没有 5.5 高。这表明该物理过程背后存在一些其他分布，可能是某种 Gamma 分布，它会在接近零的附近下降到零，但在 0-10 区间的早期达到峰值，然后有更长的尾巴。

【讨论】：

【解决方案2】：

您可以使用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)

【讨论】：