高斯混合模型与mle2 /优化

Question

我有一个mle2模型，我在这里开发只是为了演示问题。我从两个独立的高斯分布x1和x2生成值，将它们组合在一起形成x=c(x1,x2)，然后创建一个MLE，尝试将x值重新归类为属于值的左侧特定值或特定x值的右侧通过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的刻度（您可以尝试使其更加清晰;原则上，您可以将其设置为可调参数，但如果这样做，则可能会再次遇到问题，因为优化器可能希望使其成为无限...）换句话说，相当于说组件1中的所有观察结果都是肯定是，并且组件2中的所有观察结果都是肯定是，我们说观察结果是等于xsplit在任一分量中都有50/50的概率下降，随着x下降到xsplit以下，分量1中的确定性增加。具有非常大的缩放参数的逻辑函数接近先前尝试的锐分模型;一般你想让缩放参数“足够大”以得到合理的分割，并且足够小，不会遇到数字问题。 （如果你的比例太大，计算的概率会下溢/溢出到0或1，你会回到你开始的地方...）

这是我第二次或第三次尝试;我必须做相当的摆弄（边界值从0或0到1之间，并将标准偏差用对数标度拟合），但结果似乎是合理的。如果我不在逻辑（plogis）函数上使用clamp()，那么我得到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 

的结果看起来是合理的。