Obrigado colega!
E aproveitando o insejo, encontrei um artigo que aborta sobre Modelos de Mistura, o qual disponibiliza a rotina R utilizada no mesmo, entretanto, como ainda não estou bem familiarizado com o assunto, não entendi qual o objetivo da rotina. Mas de qualquer forma, a rotina é a seguinte:
library(mixtools) # mixtools package must be installed first
set.seed(123) # Ensure that results are exactly reproducible
#Generate data:
mu <- matrix(c(0, 15), 2, 3)
sigma <- matrix(c(1, 5), 2, 3)
x <- rmvnormmix(300, lambda = c(.4,.6), mu = mu, sigma = sigma)
# npEM algorithm results:
a <- npEM(x, mu0 = 2, blockid = rep(1,3), samebw = TRUE)
b <- npEM(x, mu0 = 2, blockid = rep(1,3), samebw = FALSE)
# Produce plots like in Figure 1:
s <- seq(-10, 40, len = 200)
plot(a, xlim=c(-10, 40), ylim = c(0, .16), xlab = "", breaks = 30)
lines(s, dnorm(s)*.4 + dnorm((s-15)/5)/5*.6, lwd = 2, lty = 2)
plot(b, xlim=c(-10, 40), ylim = c(0, .16), xlab = "", breaks = 30)
lines(s, dnorm(s)*.4 + dnorm((s-15)/5)/5*.6, lwd = 2, lty = 2)
# Display npEM bandwidths, minimum values of max_j p_{ij}
a$bandwidth
b$bandwidth
min(apply(a$posteriors, 1, max))
min(apply(b$posteriors, 1, max))
# spEM algorithm results:
d <- spEM(x, mu0 = 2, blockid = rep(1,3), constbw = FALSE)
d2 <- spEM(x, mu0 = 2, blockid = rep(1,3), constbw = TRUE)
# Produce plot like in Figure 2:
plot(d, xlim=c(-10, 40), ylim = c(0, .16), xlab = "", breaks = 30,
addlegend=FALSE)
plot(d2, newplot=FALSE, addlegend=FALSE, lty=2, dens.col=1)
# Display spEM bandwidths
d$bandwidth
d2$bandwidth
Caso se interesse pelo artigo, posso enviá-lo nesse mesmo e-mail.
Att.
André