Модели смеси содержат множество технических проблем (симметрия при повторной маркировке компонентов и т. Д.); если у вас нет особых потребностей, вам может быть лучше использовать одно из большого количества пакетов моделирования специального назначения, которые были написаны для R (только library("sos"); findFn("{mixture model}") или findFn("{mixture model} Gaussian")).
Однако в этом случае у вас есть более конкретная проблема, заключающаяся в том, что поверхность хорошего качества/правдоподобия параметра xsplit является «плохим» (т. Е. Производная равна нулю почти всюду). В частности, если вы считаете пару точек x1, x2 в вашем наборе данных, которые являются соседями, вероятность одинакова для любого параметра разделения между x1 и x2 (поскольку любое из этих значений разбивает набор данных на те же самые два компонента). Это означает, что поверхность правдоподобия кусочно плоская, что делает практически невозможным любой разумный оптимизатор - даже такие, как Nelder-Mead, которые явно не зависят от производных. Ваш выбор: (1) использовать какой-то стохастический оптимизатор грубой силы (например, метод = «SANN» в optim()); (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. Логистическая функция с очень большим параметром масштабирования аппроксимирует ранее опробованную модель резкого раскола; обычно вы хотите сделать параметр масштабирования «достаточно большим», чтобы получить разумный раскол и достаточно маленький, чтобы не сталкиваться с числовыми проблемами. (Если вы сделаете масштаб слишком большим, расчетные вероятности будут переполнены/переполнены до 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
Результаты выглядят разумными.
для чего это стоит, это проблема оптимизации, а не проблема 'mle2'; 'mle2' просто обертывает функцию' optim'. Модели смешивания **, как известно, трудно подстраиваются - для них разработано множество алгоритмов оптимизации специального назначения. –
, если mle2 обертывает функцию оптимизации, тогда я не понимаю, почему она объясняет, как это происходит, потому что под капотом это делает оптимизм. – CodeGuy
Было бы проще начать с использования 'nls' для соответствия функции сортировки' a1 * exp (-x^2/b1) + a2 * exp (-x^2/b2) ', а затем классифицировать данные по относительные амплитуды этих двух гауссианов? (Конечно, это не сработает, когда критерий Рэлея не очень хорошо встречен) –