r - 使用 mle2/optim 进行高斯混合建模

Question

我在这里开发了一个 mle2 模型，只是为了演示该问题。我从两个单独的高斯分布 x1 和 x2 生成值，将它们组合在一起形成 x=c(x1,x2)，然后创建尝试将 x 值重新分类为属于特定 x 值的左侧或特定 x 值的右侧的 MLE，通过xsplit 参数。

问题是找到的参数并不理想。具体来说，xsplit 始终以其起始值返回。如果我更改其起始值(例如，4 或 9)，则结果的对数似然会存在巨大差异。

这是完全可重现的示例:

set.seed(1001)
library(bbmle)
x1 = rnorm(n=100,mean=4,sd=0.8)
x2 = rnorm(n=100,mean=12,sd=0.4)
x = c(x1,x2)
hist(x,breaks=20)
ff = function(m1,m2,sd1,sd2,xsplit) {
  outs = rep(NA,length(xvals))
  for(i in seq(1,length(xvals))) {
    if(xvals[i]<=xsplit) {
      outs[i] = dnorm(xvals[i],mean=m1,sd=sd1,log=T)
    }
    else {
      outs[i] = dnorm(xvals[i],mean=m2,sd=sd2,log=T)
    }
  }
  -sum(outs)
}

# change xsplit starting value here to 9 and 4
# and realize the difference in log likelihood
# Why isn't mle finding the right value for xsplit?
mo = mle2(ff,
          start=list(m1=1,m2=2,sd1=0.1,sd2=0.1,xsplit=9), 
          data=list(xvals=x))

#print mo to see log likelihood value
mo

#plot the result
c=coef(mo)
m1=as.numeric(c[1])
m2=as.numeric(c[2])
sd1=as.numeric(c[3])
sd2=as.numeric(c[4])
xsplit=as.numeric(c[5])
leftx = x[x<xsplit]
rightx = x[x>=xsplit]
y1=dnorm(leftx,mean=m1,sd=sd1)
y2=dnorm(rightx,mean=m2,sd=sd2)
points(leftx,y1*40,pch=20,cex=1.5,col="blue")
points(rightx,y2*90,pch=20,cex=1.5,col="red")

如何修改我的 mle2 以捕获正确的参数，特别是 xsplit？

Answer 1

混合模型提出了很多技术挑战(组件重新标记下的对称性等)；除非您有非常具体的需求，否则您最好使用为 R 编写的大量专用混合建模包之一(仅 library("sos"); findFn("{mixture model}") 或 findFn("{mixture model} Gaussian") )。

但是，在这种情况下，您有一个更具体的问题，即 xsplit 的拟合优度/似然度曲面参数是“坏”的(即导数几乎在任何地方都为零)。特别是，如果您考虑一对点 x1 , x2在相邻的数据集中，x1 之间的任何分割参数的可能性完全相同。和x2 (因为任何这些值都会将数据集分成相同的两个部分)。这意味着似然曲面是分段平坦的，这使得任何明智的优化器几乎不可能实现——即使是像 Nelder-Mead 这样不明确依赖于导数的优化器。您的选择是(1)使用某种强力随机优化器(例如 optim() 中的 method="SANN")； (2)取xsplit超出您的功能并对其进行配置(即对于 xsplit 的每个可能选择，优化其他四个参数)； (3) 平滑您的分割标准(即拟合属于一个组件或另一个组件的逻辑概率)； (4) 使用上面推荐的专用混合模型拟合算法。

set.seed(1001)
library(bbmle)
x1 = rnorm(n=100,mean=4,sd=0.8)
x2 = rnorm(n=100,mean=12,sd=0.4)
x = c(x1,x2)

您的ff函数可以写得更紧凑:

## ff can be written more compactly:
ff2 <- function(m1,m2,sd1,sd2,xsplit) {
    p <- xvals<=xsplit
    -sum(dnorm(xvals,mean=ifelse(p,m1,m2),
               sd=ifelse(p,sd1,sd2),log=TRUE))
}

## ML estimation
mo <- mle2(ff2,
           start=list(m1=1,m2=2,sd1=0.1,sd2=0.1,xsplit=9), 
           data=list(xvals=x))

## refit with a different starting value for xsplit
mo2 <- update(mo,start=list(m1=1,m2=2,sd1=0.1,sd2=0.1,xsplit=4))

## not used here, but maybe handy
plotfun <- function(mo,xvals=x,sizes=c(40,90)) {
    c <- coef(mo)
    hist(xvals,col="gray")
    p <- xvals <= c["xsplit"]
    y <- with(as.list(coef(mo)),
              dnorm(xvals,mean=ifelse(p,m1,m2),
                    sd=ifelse(p,sd1,sd2))*sizes[ifelse(p,1,2)])
    points(xvals,y,pch=20,cex=1.5,col=c("blue","red")[ifelse(p,1,2)])
}

plot(slice(mo),ylim=c(-0.5,10))
plot(slice(mo2),ylim=c(-0.5,10))

我做了一点作弊，只提取了 xsplit参数:

xsplit=9 附近的似然面:

xsplit=4 附近的似然面:

另请参阅p. 243 of Bolker 2008 .

更新:平滑

正如我上面提到的，一种解决方案是使两个混合物成分之间的边界平滑或渐变，而不是尖锐。我使用了逻辑函数 plogis()中点位于xsplit和任意设置为 2 的比例(您可以尝试使其更清晰；原则上您可以将其设为可调整参数，但如果您这样做，您可能会再次遇到麻烦，因为优化器可能希望使其无限...... .) 换句话说，而不是说所有带有 x<xsplit 的观察结果肯定在组件 1 中，并且所有带有 x>xsplit 的观察值肯定在分量 2 中，我们说观测值等于 xsplit有 50/50 的概率落在任一组件中，而处于组件 1 中的确定性增加为 x跌破 xsplit 。具有非常大的缩放参数的逻辑函数近似于先前尝试的锐 split 模型；通常，您希望使缩放参数“足够大”以获得合理的分割，并且足够小以免遇到数字问题。 (如果比例太大，计算出的概率将下溢/上溢至 0 或 1，然后您将回到开始的位置...)

这是我的第二次或第三次尝试；我不得不做大量的摆弄(将值限制在远离 0 或 0 和 1 之间的范围，并在对数刻度上拟合标准差)，但结果似乎是合理的。如果我不使用clamp()在逻辑( plogis )函数上，我得到 0 或 1 概率；如果我不使用clamp() (单方面)在正态概率上，那么它们可以下溢到零——在任何一种情况下我都会得到无穷大或 NaN结果。在对数刻度上拟合标准差效果更好，因为当优化器尝试标准差为负值时不会遇到问题......

 ## bound x values between lwr and upr
 clamp <- function(x,lwr=0.001,upr=0.999) {
     pmin(upr,pmax(lwr,x))
 }

 ff3 <- function(m1,m2,logsd1,logsd2,xsplit) {
     p <- clamp(plogis(2*(xvals-xsplit)))
     -sum(log((1-p)*clamp(dnorm(xvals,m1,exp(logsd1)),upr=Inf)+
                  p*clamp(dnorm(xvals,m2,exp(logsd2)),upr=Inf)))
 }
 xvals <- x
 ff3(1,2,0.1,0.1,4)                                 
 mo3 <- mle2(ff3,
           start=list(m1=1,m2=2,logsd1=-1,logsd2=-1,xsplit=4), 
           data=list(xvals=x))
 ## Coefficients:
 ##          m1          m2      logsd1      logsd2      xsplit 
 ##  3.99915532 12.00242510 -0.09344953 -1.13971551  8.43767997 

结果看起来很合理。