R - Рассчитать MLE комбинированного распределения

У меня есть набор данных с 1000 значений, который представляет собой комбинацию двух нормальных распределений N (y1,1) и N (y2,1). Плотность выглядит следующим образом:

введите здесь описание изображения

Я хочу рассчитать долю N(y1,1) и N(y2,1) в наборе данных и два средних значения y1 и y2. Это мой текущий подход:

z <- #Dataset as vector with 1000 entries#
lik <- function(mu1, mu2, part) -sum(part*dnorm(z, mu1, 1, log=TRUE) + (1-part)*dnorm(z, mu2, 1, log=TRUE))
mle <- mle(lik, start=list(mu1=-7, mu2=5, part=0.33))

Но это дает мне следующее сообщение об ошибке:

Error in solve.default(oout$hessian) : 
    Lapack routine dgesv: system is exactly singular: U[1,1] = 0

r normal-distribution mle
person Codey    schedule 31.05.2018    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Я переопределил вероятность использования log() вместо аргумента log = TRUE.

Как ни странно, следующее работает, несмотря на предупреждения. Обратите внимание, что это предупреждения, а не ошибки.

library(stats4)

set.seed(7850)    # Make the results reproducible
z <- sample(c(rnorm(333, -7, 1), rnorm(667, 5, 1)))

plot(density(z))

lik2 <- function(mu1, mu2, part) -sum(log(part*dnorm(z, mu1, 1) + (1-part)*dnorm(z, mu2, 1)))
mle2 <- mle(lik2, start = list(mu1 = -6, mu2 = 6, part = 1/2))
#Warning messages:
#1: In log(part * dnorm(z, mu1, 1) + (1 - part) * dnorm(z, mu2, 1)) :
#  NaNs produced
#2: In log(part * dnorm(z, mu1, 1) + (1 - part) * dnorm(z, mu2, 1)) :
#  NaNs produced
#3: In log(part * dnorm(z, mu1, 1) + (1 - part) * dnorm(z, mu2, 1)) :
#  NaNs produced
#4: In log(part * dnorm(z, mu1, 1) + (1 - part) * dnorm(z, mu2, 1)) :
#  NaNs produced

mle2
#
#Call:
#mle(minuslogl = lik2, start = list(mu1 = -6, mu2 = 6, part = 1/2))
#
#Coefficients:
#       mu1        mu2       part 
#-7.1091780  4.9377339  0.3330038
person Rui Barradas    schedule 31.05.2018
comment
Спасибо, работает идеально. Странно, что это не работает с log = TRUE - person Codey; 31.05.2018
comment
@Codey Обратите внимание, что log(a*dnorm + b*dnorm) не равно a*log(dnorm) + b*log(dnorm). - person Rui Barradas; 31.05.2018
comment
О, точно, как-то не подумал об этом. Так что математическое log(a*dnorm + b*dnorm) в любом случае является правильным решением. Спасибо - person Codey; 31.05.2018