1

我正在尝试使用 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 视为“赢家”（因为我产生了一个单一的正态分布）。但是，如您所见，我得到 1 的无限值，而 BIC、ICL 和 Loglik 的 2 和 3 具有相同的值。这违背了正态分布和更高（或相同）的概率来处理双或多模式分布。由于我一开始使用的是正态分布，因此我希望看到 1 的概率最高，并且至少在 2 和 3 之间存在一些差异。最让我困惑的是在某些测试中 2 和 3 的相同值。

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

4

1 回答 1

0

我对mixtools包裹了解不多。我确实对您所做的进行了一些尝试，但我没有得出与您相同的结论。

当我拟合一个双分量多项式混合模型（这就是你正在做的multmixmodel.sel）时，第二个分量是不存在的；后验概率几乎为零。

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)
于 2018-02-27T12:22:25.320 回答