按一个因子将 MASS::fitdistr 应用于多个数据

Question

我的问题以粗体结尾。

我知道如何将 beta 分布拟合到某些数据。例如：

library(Lahman)
library(dplyr)

# clean up the data and calculate batting averages by playerID
batting_by_decade <- Batting %>%
  filter(AB > 0) %>%
  group_by(playerID, Decade = round(yearID - 5, -1)) %>%
  summarize(H = sum(H), AB = sum(AB)) %>%
  ungroup() %>%
  filter(AB > 500) %>%
  mutate(average = H / AB)

# fit the beta distribution
library(MASS)
m <- MASS::fitdistr(batting_by_decade$average, dbeta,
                    start = list(shape1 = 1, shape2 = 10))

alpha0 <- m$estimate[1]
beta0 <- m$estimate[2]

# plot the histogram of data and the beta distribution
ggplot(career_filtered) +
  geom_histogram(aes(average, y = ..density..), binwidth = .005) +
  stat_function(fun = function(x) dbeta(x, alpha0, beta0), color = "red",
                size = 1) +
  xlab("Batting average")

产生：

现在我想为数据的每个 batting_by_decade$Decade 列计算不同的 beta 参数 alpha0 和 beta0，所以我最终得到 15 个参数集和 15 个 beta 分布，我可以拟合这个击球平均值的 ggplot由十年分面：

batting_by_decade %>% 
  ggplot() +
  geom_histogram(aes(x=average)) +
  facet_wrap(~ Decade)

我可以通过过滤每个十年来硬编码，并将该十年的数据传递给fidistr 函数，在所有十年中重复此操作，但是有没有一种方法可以快速计算每十年的所有 beta 参数并且可重现，也许使用其中一种应用功能？

Answer 1

您可以将summarise 与两个自定义函数一起使用：

getAlphaEstimate = function(x) {MASS::fitdistr(x, dbeta,start = list(shape1 = 1, shape2 = 10))$estimate[1]}

getBetaEstimate = function(x) {MASS::fitdistr(x, dbeta,start = list(shape1 = 1, shape2 = 10))$estimate[2]}

batting_by_decade %>%
  group_by(Decade) %>%
  summarise(alpha = getAlphaEstimate(average),
         beta = getBetaEstimate(average)) -> decadeParameters

但是，根据 Hadley 的帖子，您将无法使用 stat_summary 绘制它：https://stackoverflow.com/a/1379074/3124909

Answer 2

这是一个示例，说明您如何从生成虚拟数据一直到绘图。

temp.df <- data_frame(yr = 10*187:190,
                      al = rnorm(length(yr), mean = 4, sd = 2),
                      be = rnorm(length(yr), mean = 10, sd = 2)) %>% 
  group_by(yr, al, be) %>% 
  do(data_frame(dats = rbeta(100, .$al, .$be)))

首先我制定了一些四年的规模参数，按每个组合分组，然后使用do 创建一个包含来自每个分布的 100 个样本的数据框。除了知道“真实”参数之外，这个数据框应该看起来很像您的原始数据：具有相关年份的样本向量。

---

temp.ests <- temp.df %>% 
  group_by(yr, al, be) %>% 
  summarise(ests = list(MASS::fitdistr(dats, dbeta, start = list(shape1 = 1, shape2 = 1))$estimate)) %>% 
  unnest %>% 
  mutate(param = rep(letters[1:2], length(ests)/2)) %>% 
  spread(key = param, value = ests)

这是您的大部分问题，以您解决问题的方式解决了很多问题。如果你逐行遍历这个 sn-p，你会看到你有一个数据框，其中有一列类型为list，每行都包含<dbl [2]>。当您unnest() 时，它会将这两个数字分成单独的行，因此我们通过添加一列“a，b，a，b，...”和 spread 将它们分开以获得两列来识别它们每年一排。在这里，您还可以看到fitdistr 与我们抽样的真实人口的匹配程度，查看a 与al 和b 与be。

---

temp.curves <- temp.ests %>% 
  group_by(yr, al, be, a, b) %>% 
  do(data_frame(prop = 1:99/100,
                trueden = dbeta(prop, .$al, .$be),
                estden = dbeta(prop, .$a, .$b)))

现在我们将这个过程彻底翻转，以生成数据来绘制曲线。对于每一行，我们使用do 制作一个包含一系列值prop 的数据框，并计算真实总体参数和我们估计的样本参数在每个值处的β 密度。

---

ggplot() +
  geom_histogram(data = temp.df, aes(dats, y = ..density..), colour = "black", fill = "white") +
  geom_line(data = temp.curves, aes(prop, trueden, color = "population"), size = 1) +
  geom_line(data = temp.curves, aes(prop, estden, color = "sample"), size = 1) +
  geom_text(data = temp.ests, 
            aes(1, 2, label = paste("hat(alpha)==", round(a, 2))), 
            parse = T, hjust = 1) +
  geom_text(data = temp.ests, 
            aes(1, 1, label = paste("hat(beta)==", round(b, 2))), 
            parse = T, hjust = 1) +
  facet_wrap(~yr)

最后，我们将它们放在一起，绘制了样本数据的直方图。然后是我们的曲线数据中的一条线作为真实密度。然后是我们估计密度的曲线数据中的一条线。然后我们的参数估计数据中的一些标签来显示样本参数，以及按年份显示的方面。

Answer 3

这是一个应用解决方案，但我更喜欢@CMichael 的 dplyr 解决方案。

calc_beta <- function(decade){
  dummy <- batting_by_decade %>% 
    dplyr::filter(Decade == decade) %>% 
    dplyr::select(average)

  m <- fitdistr(dummy$average, dbeta, start = list(shape1 = 1, shape2 = 10))

  alpha0 <- m$estimate[1]
  beta0 <- m$estimate[2]

  return(c(alpha0,beta0))
}

decade <- seq(1870, 2010, by =10)
params <- sapply(decade, calc_beta)
colnames(params) <- decade

回复：@CMichael 关于避免双重fitdistr 的评论，我们可以将函数重写为getAlphaBeta。

getAlphaBeta = function(x) {MASS::fitdistr(x, dbeta,start = list(shape1 = 1, shape2 = 10))$estimate}

batting_by_decade %>%
  group_by(Decade) %>%
  summarise(params = list(getAlphaBeta(average))) -> decadeParameters

decadeParameters$params[1] # it works!

现在我们只需要以一种不错的方式取消列出第二列....