通過一個因子將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") 

其中產量：

現在我想計算不同的測試參數alpha0和beta0對數據的每個batting_by_decade$Decade列，這樣我結束了15個參數集，和15個分佈，我可以適應這個ggplot打擊平均面Decade：

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

我可以通過過濾每一個十年，並通過數據的十年的身價進入fidistr功能，重複此爲所有幾十年，但硬編碼，這是有快速計算每十年所有測試參數的方法並可重複，也許與其中一個應用功能？

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 

但是，您將無法根據哈德利的帖子在這裏與stat_summary繪製它：https://stackoverflow.com/a/1379074/3124909

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個樣本的數據幀。除了知道「真實」參數之外，這個數據框應該看起來很像您的原始數據：具有相關年份的樣本矢量。

---

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 vs al和b vs 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) 

最後，我們把它放在一起，密謀我們的樣本數據的直方圖。然後從我們的曲線數據中獲得真實密度的一條線。然後從我們的曲線數據中獲得一條線，用於估算密度。然後根據我們的參數估計數據中的一些標籤來顯示樣本參數，以及按年份顯示的方面。