Beta 分布的混合：完整示例

机器算法验证 r 贝塔回归 百分比
2022-03-16 11:00:10

我有一个测量比例的向量并且想要

  1. 测试这些比例是否遵循双峰分布
  2. 表征两个基础分布
  3. 确定每个数据点属于哪个分布
  4. 总体上直方图上的密度分布

我一直在使用 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

在此处输入图像描述

1个回答

问题中的所有项目都没有完全正确解决。我会推荐以下。

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

结果图如下所示。

在此处输入图像描述

其它你可能感兴趣的问题
上一篇退化的单变量高斯 下一篇在 OLS 回归中使用非平稳时间序列数据