正态分布对多模态的概率高于单峰态
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)
我正在尝试使用 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)