[Depmix-commits] r533 - in pkg/depmixS4: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 22 19:21:34 CEST 2012
Author: maarten
Date: 2012-06-22 19:21:33 +0200 (Fri, 22 Jun 2012)
New Revision: 533
Modified:
pkg/depmixS4/R/EM.R
pkg/depmixS4/R/depmixfit.R
pkg/depmixS4/R/em.control.R
pkg/depmixS4/man/em.control.Rd
Log:
- implementation of classification likelihood in EM complete (still needs thorough testing)
Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R 2012-06-21 17:55:49 UTC (rev 532)
+++ pkg/depmixS4/R/EM.R 2012-06-22 17:21:33 UTC (rev 533)
@@ -59,7 +59,6 @@
if(clsf == "hard") {
gamma <- t(apply(gamma,1,ind.max))
}
-
LL <- -1e10
for(i in 1:ns) {
@@ -71,7 +70,13 @@
}
# initial expectation
- fbo <- fb(init=object at init,matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
+ if(clsf == "hard") {
+ fbo <- list()
+ vstate <- apply(gamma,1,which.max)
+ fbo$logLike <- sum(log((apply(object at dens,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+ } else {
+ fbo <- fb(init=object at init,matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
+ }
LL <- fbo$logLike
if(is.nan(LL)) stop("Cannot find suitable starting values; please provide them.")
@@ -79,10 +84,18 @@
} else {
# initial expectation
B <- apply(object at dens,c(1,3),prod)
- gamma <- object at init*B
- LL <- sum(log(rowSums(gamma)))
+
+ if(clsf == "hard") {
+ gamma <- t(apply(object at init*B,1,ind.max))
+ LL <- sum(log(rowSums(gamma*B)))
+ } else {
+ gamma <- object at init*B
+ LL <- sum(log(rowSums(gamma)))
+ gamma <- gamma/rowSums(gamma)
+ }
+
if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
- gamma <- gamma/rowSums(gamma)
+
}
LL.old <- LL + 1
@@ -93,9 +106,6 @@
# should become object at prior <- fit(object at prior)
- if(clsf == "hard") {
- gamma <- t(apply(gamma,1,ind.max))
- }
object at prior@y <- gamma[bt,,drop=FALSE]
object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
object at init <- dens(object at prior)
@@ -110,12 +120,20 @@
# expectation
B <- apply(object at dens,c(1,3),prod)
- gamma <- object at init*B
- LL <- sum(log(rowSums(gamma)))
+ if(clsf == "hard") {
+ gamma <- t(apply(object at init*B,1,ind.max))
+ LL <- sum(log(rowSums(gamma*B)))
+ } else {
+ gamma <- object at init*B
+ LL <- sum(log(rowSums(gamma)))
+ # normalize
+ gamma <- gamma/rowSums(gamma)
+ }
+ #gamma <- object at init*B
+
- # normalize
- gamma <- gamma/rowSums(gamma)
+
# print stuff
if(verbose&((j%%5)==0)) {
cat("iteration",j,"logLik:",LL,"\n")
@@ -139,10 +157,17 @@
class(object) <- "mix.fitted"
if(converge) {
- object at message <- switch(crit,
- relative = "Log likelihood converged to within tol. (relative change)",
- absolute = "Log likelihood converged to within tol. (absolute change)"
- )
+ if(clsf == "hard") {
+ object at message <- switch(crit,
+ relative = "Log classification likelihood converged to within tol. (relative change)",
+ absolute = "Log classification likelihood converged to within tol. (absolute change)"
+ )
+ } else {
+ object at message <- switch(crit,
+ relative = "Log likelihood converged to within tol. (relative change)",
+ absolute = "Log likelihood converged to within tol. (absolute change)"
+ )
+ }
} else object at message <- "'maxit' iterations reached in EM without convergence."
# no constraints in EM, except for the standard constraints ...
@@ -193,14 +218,32 @@
}
# initial expectation
- fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ if(clsf == "hard") {
+ fbo <- list()
+ vstate <- viterbi(object)[,1]
+ fbo$gamma <- as.matrix(model.matrix(~ factor(vstate) - 1))
+ fbo$xi <- array(0,dim=c(sum(ntimes),ns,ns))
+ fbo$xi[cbind(1:(sum(ntimes)- 1),vstate[-1],vstate[-length(vstate)])] <- 1
+ fbo$logLike <- sum(log((apply(object at dens,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+ } else {
+ fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ }
LL <- fbo$logLike
if(is.nan(LL)) stop("Cannot find suitable starting values; please provide them.")
} else {
# initial expectation
- fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ if(clsf == "hard") {
+ fbo <- list()
+ vstate <- viterbi(object)[,1]
+ fbo$gamma <- as.matrix(model.matrix(~ factor(vstate) - 1))
+ fbo$xi <- array(0,dim=c(sum(ntimes),ns,ns))
+ fbo$xi[cbind(1:(sum(ntimes)- 1),vstate[-1],vstate[-length(vstate)])] <- 1
+ fbo$logLike <- sum(log((apply(object at dens,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+ } else {
+ fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ }
LL <- fbo$logLike
if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
}
@@ -212,16 +255,7 @@
# maximization
# should become object at prior <- fit(object at prior, gamma)
-
- if(clsf == "hard") {
- vstate <- viterbi(object)[,1]
- fbo$gamma <- as.matrix(model.matrix(~ factor(vstate) - 1))
- # TODO: compute fbo$xi
- fbo$xi <- array(0,dim=dim(fbo$xi))
- fbo$xi[cbind(1:(dim(fbo$xi)[1] - 1),vstate[-1],vstate[-length(vstate)])] <- 1
- # TODO: check likelihood
- }
-
+
object at prior@y <- fbo$gamma[bt,,drop=FALSE]
object at prior <- fit(object at prior, w=NULL, ntimes=NULL)
object at init <- dens(object at prior)
@@ -256,10 +290,18 @@
}
}
- # expectation
- fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
- LL <- fbo$logLike
-
+ if(clsf == "hard") {
+ vstate <- viterbi(object)[,1]
+ fbo$gamma <- as.matrix(model.matrix(~ factor(vstate) - 1))
+ fbo$xi <- array(0,dim=c(sum(ntimes),ns,ns))
+ fbo$xi[cbind(1:(sum(ntimes)- 1),vstate[-1],vstate[-length(vstate)])] <- 1
+ fbo$logLike <- sum(log((apply(object at dens,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+ } else {
+ # expectation
+ fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ }
+
+ LL <- fbo$logLike
if(verbose&((j%%5)==0)) cat("iteration",j,"logLik:",LL,"\n")
if( (LL >= LL.old)) {
Modified: pkg/depmixS4/R/depmixfit.R
===================================================================
--- pkg/depmixS4/R/depmixfit.R 2012-06-21 17:55:49 UTC (rev 532)
+++ pkg/depmixS4/R/depmixfit.R 2012-06-22 17:21:33 UTC (rev 533)
@@ -29,7 +29,7 @@
if(method=="EM") {
if(!(emcontrol$crit %in% c("absolute","relative"))) stop("'crit' argument to em.control not recognized")
- object <- em(object,maxit=emcontrol$maxit,tol=emcontrol$tol,crit=emcontrol$crit,random.start=emcontrol$random.start,verbose=verbose,...)
+ object <- em(object,maxit=emcontrol$maxit,tol=emcontrol$tol,crit=emcontrol$crit,random.start=emcontrol$random.start,classification=emcontrol$classification,verbose=verbose,...)
}
if(!(method %in% c("EM","donlp","rsolnp"))) stop("'method' argument invalid; should be one of 'EM', 'rsolnp', 'donlp'.")
Modified: pkg/depmixS4/R/em.control.R
===================================================================
--- pkg/depmixS4/R/em.control.R 2012-06-21 17:55:49 UTC (rev 532)
+++ pkg/depmixS4/R/em.control.R 2012-06-22 17:21:33 UTC (rev 533)
@@ -1,4 +1,5 @@
em.control <-
-function(maxit=500,tol=1e-8,crit="relative",random.start=TRUE,classification=c("soft","hard")) {
- classification <- match.arg(classification)
+function(maxit=500,tol=1e-8,crit=c("relative","absolute"),random.start=TRUE,classification=c("soft","hard")) {
+ crit <- match.arg(crit)
+ classification <- match.arg(classification)
return(list(maxit=maxit,tol=tol,crit=crit,random.start=random.start,classification=classification))}
Modified: pkg/depmixS4/man/em.control.Rd
===================================================================
--- pkg/depmixS4/man/em.control.Rd 2012-06-21 17:55:49 UTC (rev 532)
+++ pkg/depmixS4/man/em.control.Rd 2012-06-22 17:21:33 UTC (rev 533)
@@ -8,7 +8,7 @@
\usage{
- em.control(maxit = 500, tol = 1e-08, crit = "relative", random.start = TRUE, classification = c("soft","hard"))
+ em.control(maxit = 500, tol = 1e-08, crit = c("relative","absolute"), random.start = TRUE, classification = c("soft","hard"))
}
More information about the depmix-commits
mailing list