按一个因子将 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 分布可以适合这个由 Decade 分面的击球平均值的 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 的 post，您将无法使用 stat_summary 绘制它：

Answer 2

这是一个应用解决方案，但我更喜欢@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!

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

Answer 3

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

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 个样本。除了了解 "true" 参数之外，此数据框应该看起来很像您的原始数据：具有相关年份的样本向量。

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)

这是你这里的大部分问题，按照你解决的方式解决了很多问题。如果逐行执行此代码段，您将看到一个数据框，其中有一列类型为 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 的数据框，并计算每个值的真实总体参数和我们估计的样本参数的 beta 密度。

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)

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

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

Apply MASS::fitdistr to multiple data by a factor

r

apply

ggplot2

beta-distribution

dplyr