2

我的问题最后以粗体显示。

我知道如何将 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")

产生：

在此处输入图像描述

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

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

在此处输入图像描述

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

4

3 回答 3

2

您可以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

但是，您将无法stat_summary根据 Hadley 在此处的帖子进行绘制：https ://stackoverflow.com/a/1379074/3124909

于 2017-08-12T20:06:50.200 回答
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)

这是您的大部分问题，非常解决了您解决问题的方式。如果您逐行浏览此代码段，您会看到您有一个数据框，其中有一列 type list，包含<dbl [2]>在每一行中。当您unnest()将这两个数字拆分为单独的行时，然后我们通过添加一列“a，b，a，b，...”来识别它们，然后将spread它们分开以获得每年一行的两列。在这里，您还可以通过查看vs和vs来了解fitdistr我们从中抽样的真实总体的匹配程度。aalbbe

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)

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

在此处输入图像描述

于 2017-08-13T03:31:14.170 回答
1

这是一个应用解决方案，但我更喜欢@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 关于避免使用 double 的评论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!

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

于 2017-08-12T20:26:55.740 回答