4

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

问题是找到的参数并不理想。具体来说，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？

4

1 回答 1

8

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

但是，在这种情况下，您有一个更具体的问题，即参数的拟合优度/似然面xsplit是“坏的”（即导数几乎处处为零）。特别是，如果您考虑数据集中的一对相邻点x1，则和x2之间的任何拆分参数的可能性完全相同（因为这些值中的任何一个都将数据集拆分为相同的两个组件）。这意味着似然面是分段平坦的，这使得任何明智的优化器几乎都不可能——即使是像 Nelder-Mead 这样不明确依赖导数的优化器。您的选择是（1）使用某种强力随机优化器（例如 optim() 中的 method="SANN"）；(2) 取x1x2xsplit超出您的功能和配置文件（即对于每个可能的选择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=9

周围的似然面xsplit=4：

xsplit=4

另见p. Bolker 2008 年第 243 页。

更新：平滑

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

结果看起来很合理。

于 2014-02-07T17:35:49.893 回答