regmixEM.loc=function (y, x, lambda = NULL, beta = NULL, sigma = NULL, k = 2,
addintercept = TRUE, kern.l = c("Gaussian", "Beta", "Triangle",
"Cosinus", "Optcosinus"), epsilon = 1e-08,
maxit = 10000, kernl.g = 0, kernl.h = 1, verb=FALSE)
{
diff <- 1
iter <- 0
x.1 <- cbind(1, x)
n <- length(y)
kern.l <- match.arg(kern.l)
out.EM <- regmixEM.lambda(y, x, lambda = lambda, beta = beta,
sigma = sigma, k = k, epsilon = epsilon, maxit = maxit)
ll = out.EM$loglik
restarts <- 0
while (diff > epsilon && iter < maxit) {
old.l.x = out.EM$lambda
old.loglik = out.EM$loglik
l.x = matrix(nrow = n, ncol = (k - 1))
for (i in 1:(k - 1)) {
l.x[, i] <- apply(matrix(x, ncol = 1), 1, lambda,
z = out.EM$post[, i], xi = x, kernel = kern.l,
g = kernl.g, h = kernl.h)
}
l.x <- cbind(l.x, 1 - apply(l.x, 1, sum))
out.EM.loc <- regmixEM.lambda(y, x, beta = out.EM$beta,
sigma = out.EM$sigma, lambda = l.x, k = k)
loglik.loc<-out.EM.loc$loglik
out.EM.old <- regmixEM.lambda(y, x, beta = out.EM$beta,
sigma = out.EM$sigma, lambda = old.l.x, k = k)
loglik.old<-out.EM.old$loglik
if(loglik.loc>old.loglik){
out.EM <- out.EM.loc
} else out.EM <- out.EM.old
loglik.chosen <- out.EM$loglik
ll <- c(ll, loglik.chosen)
diff <- loglik.chosen - old.loglik
if(diff < 0){
cat("Generating new starting values...","\n")
out.EM <- regmixEM.lambda(y, x, lambda = lambda, beta = beta,
sigma = sigma, k = k, epsilon = epsilon, maxit = maxit)
restarts <- restarts + 1
if (restarts > 15)
stop("Too many tries!")
iter<-0
diff<-1
} else{
iter <- iter + 1
if (verb) {
cat("iteration=", iter, "diff=", diff, "log-likelihood",
loglik.chosen, "\n")
}
}
}
if (iter == maxit) {
cat("WARNING! NOT CONVERGENT!", "\n")
}
cat("number of overall iterations=", iter, "\n")
a=list(x=x, y=y, lambda.x = out.EM$lambda, beta = out.EM$beta, sigma = out.EM$sigma,
loglik = loglik.chosen, posterior = out.EM$post, all.loglik=ll, restarts=restarts,
ft="regmixEM.loc")
class(a) = "mixEM"
a
}