高斯混合模型與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 

的結果看起來是合理的。