正态分布对多模态的概率高于单峰态

Normal distribution has higher probability for multimodality than unimodality

我正在尝试使用 mixtool 包分析混合模型,换句话说,我想分析我的数据是单峰分布、双峰分布还是多峰分布。

为简单起见,这里举个例子:

library(mixtools)
#creating an aritifical normal distribution
mydata <- rnorm(1000, 1750, 60)

#defining the cuts and preparing it for calculations
cutp <- seq(1600, 2300, by=25)
mult <- makemultdata(mydata, cuts = cutp)
comp <- multmixmodel.sel(mult, comps = 1:3, epsilon = 0.01)

#plotting the data (in this case 2 subpopulations)
mixmdl = normalmixEM(mydata, k=2, maxit=50000)
plot(mixmdl,which=2)
lines(density(mydata), lty=2, lwd=2)

现在 'comp' 的结果是:

          1         2          3 Winner
AIC    -Inf -94.04097 -124.04097      2
BIC    -Inf -35.04097  -35.04097      2
CAIC   -Inf -64.54097  -79.54097      2
ICL    -Inf -35.04097  -35.04097      2
Loglik -Inf -35.04097  -35.04097      2

根据我对这种执行的非常有限的理解,我希望将 1 视为 'winner'(因为我生成了单一正态分布)。 但是,如您所见,我得到 1 的无穷大值,以及 2 和 3 的 BIC、ICL 和 Loglik 的相同值。这说明正态分布和更高(或相同)的概率处理双或多模态分布。由于我一开始使用的是正态分布,因此我希望看到 1 的概率最高,并且 2 和 3 之间至少存在一些差异。最让我困惑的是某些测试中 2 和 3 的相同值。

所以我的问题是为什么我的方法无法将分布识别为高斯分布,而是将其分类为双峰/多峰分布?

我对 mixtools 软件包了解不多。我确实对你所做的做了一点尝试,但我没有得出和你一样的结论。

当我拟合一个 two-component 多项式混合模型时(这就是你用 multmixmodel.sel 所做的),第二个分量是 non-existent;后验概率几乎为零。

set.seed(1)
mydata <- rnorm(1000, 1750, 60)

cutp <- seq(min(mydata), max(mydata), by=25)
mult <- makemultdata(mydata, cuts = cutp)

multmod2 <- multmixEM(mult, k=2)

multmod2 $posterior
# comp.1        comp.2
# [1,]      1 1.980052e-226

当我将混合模型拟合到原始数据时,每次都选择单个组件。

library(mclust)
fit <- Mclust(mydata)
fit
#'Mclust' model object:
#  best model: univariate normal (X) with 1 components


library(EMMIX)
# Available from 
#https://people.smp.uq.edu.au/GeoffMcLachlan/mix_soft/EMMIX_R/

fit_1 <- EMMIX(mydata, g=1)
fit_2 <- EMMIX(mydata, g=2)

c(fit_1$bic, fit_2$bic)
# [1] 11108.02 11128.67
#(BIC selects the one component model)