MASS :: fitdistrを複数のデータに係数で適用する

Question

私の質問は末尾がで太字になっています。

私はいくつかのデータにベータ版を適合させる方法を知っています。

得

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") 

：例えば

を今、私は15枚のパラメータセットで終わるように、データの各batting_by_decade$Decade列に異なるベータパラメータにalpha0とbeta0を計算すると、 15歳のベータ版で、10年前のこの打球の平均値に合わせることができます：

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

私は、各十年のためのフィルタリング、およびfidistr関数にデータの十年の価値を渡して、すべての数十年のためにこれを繰り返しますが、によってハードコードこれはすぐに十年ごとにすべてのベータパラメータを算出する方法がありますすることができますおそらく適用関数の1つを使って再現可能になるでしょうか？

Answer 1

あなたはこのために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でそれをプロットすることはできません。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 

日時：ダブルfitdistrを回避についてCMichaelさんのコメント@、我々は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))) 

まず私は、各組み合わせごとにグループ化され、4年間、いくつかのスケールパラメータを作り、その後、各分布からの100個のサンプルを有するデータフレームを作成するためにdoを使用します。 「真の」パラメータを知ることとは別に、このデータフレームは元のデータとよく似ているはずです。

---

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) 

これは非常にあなたがそれを解決する方法を解決し、ここにあなたの質問の大半です。このスニペットを1行ずつ進めば、の列のデータフレームがあり、各行に<dbl [2]>が含まれていることがわかります。 unnest()は、これらの2つの数値を別々の行に分割するので、「a、b、a、b、...」になる列を追加してそれらを識別し、spreadを戻してそれぞれを1行ずつ2列にします年。ここでは、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) 

最後に、我々は、我々のサンプルデータのヒストグラムをプロットする、一緒にそれを置きます。その後、真の密度についての我々の曲線データからの線。次に、推定された密度についての曲線データからの線。その後、私たちのパラメータからのいくつかのラベルは、サンプルパラメータと1年ごとのファセットを示すデータを見積もります。