Beta 分佈的混合：完整示例
我有一個測量比例的向量並且想要
- 測試這些比例是否遵循雙峰分佈
- 表徵兩個基礎分佈
- 確定每個數據點屬於哪個分佈
- 總體上直方圖上的密度分佈
我一直在使用 betareg 包中的“betamix”，但還沒有弄清楚第 2 步和第 3 步。我在 stackoverflow 中搜索了解決方案，但沒有找到明確的答案。
這是我當前的代碼：
# parameters of distribution #1 alpha1 <- 10 beta1 <- 30 # parameters of distribution #2 alpha2 <- 30 beta2 <- 10 # Generate bimodal data set.seed(0) d <- data.frame(y = c(rbeta(100, alpha1, beta1), rbeta(50, alpha2, beta2))) # correction recommended in cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf, cf first paragraph in section 2, page 3): n <- length(d$y) d$yc <- (d$y* (n-1)+0.5)/n # histogram par(mfrow=c(1,2)) hist(d$yc, 50) # fitting mixtures of beta distributions uni.modal <- betamix(yc ~ 1 | 1, data = d, k = 1) bi.modal <- betamix(yc ~ 1 | 1, data = d, k = 2) # 1) test for bimodality: lrtest(uni.modal, bi.modal) # 3) determine for each data point which distribution it belongs d$group <- (posterior(bi.modal)[,1] <= posterior(bi.modal)[,2])+1 plot(d$y, col=d$group) # 2) characterize the two underlying distributions # betamix uses a different parametrization than dbeta (cf cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf). # instead of alpha and beta, betareg parametrizes the beta distribution using mu and phi, where # mu=alpha/(alpha+beta) # phi= alpha + beta # ERROR: converting mu and phi to alpha and beta mu1 <- coef(bi.modal)[1,1] phi1 <- coef(bi.modal)[1,2] (a1 <- mu1*phi1) # output: 4.61028 (b1 <- (1-mu1)*phi1) # ouput:-0.7240308 # 4) overaly the density distributions on the histogram
問題中的所有項目都沒有完全正確解決。我會推薦以下。
(0) 觀測值“y”不需要更正，因為它們已經在 0 和 1 之間。應用更正應該不會產生問題，但也沒有必要。
(1) 不能通過似然比 (LR) 檢驗來回答。通常在混合模型中，成分數量的選擇不能基於 LR 檢驗，因為它的正則性假設沒有得到滿足。相反，經常使用信息標準，並且 betamix() 所基於的“flexmix”提供 AIC、BIC 和 ICL。因此，您可以通過以下方式在 1、2、3 個集群中選擇最佳 BIC 解決方案
library("flexmix") set.seed(0) m <- betamix(y ~ 1 | 1, data = d, k = 1:3)(2) betamix() 中的參數不是直接的 mu 和 phi，而是兩個參數都使用了額外的鏈接函數。默認值分別為 logit 和 log。這可確保參數分別在其有效範圍 (0, 1) 和 (0, inf) 內。可以修改兩個組件中的模型，以便更輕鬆地訪問鏈接和反向鏈接等。但是，在這里手動應用反向鏈接可能是最簡單的：
mu <- plogis(coef(m)[,1]) phi <- exp(coef(m)[,2])這表明均值非常不同（0.25 和 0.77），而精度卻非常相似（49.4 和 47.8）。然後我們可以轉換回 alpha 和 beta，得到 12.4、37.0 和 36.7、11.1，這與模擬中的原始參數相當接近：
a <- mu * phi b <- (1 - mu) * phi(3) 可以使用 clusters() 函數提取簇。這只是選擇具有最高後驗（）概率的組件。在這種情況下，terior() 非常明確，即接近於零或接近於 1。
cl <- clusters(m)(4) 當用直方圖可視化數據時，可以分別可視化兩個分量，即每個分量都有自己的密度函數。或者可以繪製一個具有相應關節密度的關節直方圖。不同之處在於後者需要考慮不同的集群大小：這裡的先驗權重約為 1/3 和 2/3。可以像這樣繪製單獨的直方圖：
## separate histograms for both clusters hist(subset(d, cl == 1)$y, breaks = 0:25/25, freq = FALSE, col = hcl(0, 50, 80), main = "", xlab = "y", ylim = c(0, 9)) hist(subset(d, cl == 2)$y, breaks = 0:25/25, freq = FALSE, col = hcl(240, 50, 80), main = "", xlab = "y", ylim = c(0, 9), add = TRUE) ## lines for fitted densities ys <- seq(0, 1, by = 0.01) lines(ys, dbeta(ys, shape1 = a[1], shape2 = b[1]), col = hcl(0, 80, 50), lwd = 2) lines(ys, dbeta(ys, shape1 = a[2], shape2 = b[2]), col = hcl(240, 80, 50), lwd = 2) ## lines for corresponding means abline(v = mu[1], col = hcl(0, 80, 50), lty = 2, lwd = 2) abline(v = mu[2], col = hcl(240, 80, 50), lty = 2, lwd = 2)和聯合直方圖：
p <- prior(m$flexmix) hist(d$y, breaks = 0:25/25, freq = FALSE, main = "", xlab = "y", ylim = c(0, 4.5)) lines(ys, p[1] * dbeta(ys, shape1 = a[1], shape2 = b[1]) + p[2] * dbeta(ys, shape1 = a[2], shape2 = b[2]), lwd = 2)結果圖如下所示。
