0

我有以下数据集：

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

4

2 回答 2

1

您正在处理一个分类变量“间隔”，它根据假定的基础连续变量创建一个离散的计数观察，您从中获取断点。那种凌乱的数据情况。从技术上讲，您有区间删失数据。但是，如果您有指数分布作为假设，那么您计算的那些“平均值”实际上是中点，但它们不会被期望是指数分布变量的平均值。有关我对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 高。这表明其他一些分布是这个物理过程的基础，也许是某种伽玛分布，它会在接近零的时候降到零，但在 0-10 区间的早期达到峰值，然后有更长的尾巴。

于 2020-02-27T02:06:40.767 回答
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)
于 2020-02-27T00:51:38.660 回答