[Depmix-commits] r518 - in pkg/depmixS4: R man

[noreply at r-forge.r-project.org]

Author: maarten
Date: 2012-06-14 20:46:05 +0200 (Thu, 14 Jun 2012)
New Revision: 518

Modified:
   pkg/depmixS4/R/EM.R
   pkg/depmixS4/R/em.control.R
   pkg/depmixS4/man/em.control.Rd
Log:
- started with classification likelihood (done for mix models, not for depmix models)
- still needs to be tested!

Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R	2012-06-12 14:31:19 UTC (rev 517)
+++ pkg/depmixS4/R/EM.R	2012-06-14 18:46:05 UTC (rev 518)
@@ -10,7 +10,19 @@
     x/as.vector(sm)
 }
 
+which.is.max <- function(x) {
+    # taken from MASS
+    y <- seq_along(x)[x == max(x)]
+    if(length(y) > 1L) sample(y, 1L) else y
+}
 
+
+ind.max <- function(x) {
+    out <- rep(0,length(x))
+    out[which.is.max(x)] <- 1
+    out
+}
+
 em <- function(object,...) {
 	if(!is(object,"mix")) stop("object is not of class '(dep)mix'")
 	call <- match.call()
@@ -24,8 +36,10 @@
 }
 
 # em for lca and mixture models
-em.mix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,...) {
+em.mix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,classification=c("soft","hard"),...) {
 	
+	clsf <- match.arg(classification)
+	
 	if(!is(object,"mix")) stop("object is not of class 'mix'")
 		
 	ns <- nstates(object)
@@ -40,15 +54,17 @@
 	if(random.start) {
 				
 		nr <- sum(ntimes(object))
-		#gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nrow=nr,ncol=ns)
-		#gamma <- matrix(rbeta(nr*ns,alpha=.1,beta=.1),nrow=nr,ncol=ns)
 		gamma <- rdirichlet(nr,alpha=rep(.01,ns))
-		#gamma <- gamma/rowSums(gamma)
+		
+		if(clsf == "hard") {
+		    gamma <- t(apply(gamma,1,ind.max))
+		}
+
 		LL <- -1e10
 		
 		for(i in 1:ns) {
 			for(k in 1:nresp(object)) {
-				object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=gamma[,i])
+			    object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=gamma[,i])
 				# update dens slot of the model
 				object at dens[,k,i] <- dens(object at response[[i]][[k]])
 			}
@@ -76,6 +92,10 @@
 		# maximization
 		
 		# 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)
@@ -157,6 +177,9 @@
 		#gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nrow=nr,ncol=ns)
 		#gamma <- gamma/rowSums(gamma)
 		gamma <- rdirichlet(nr,alpha=rep(.01,ns))
+		if(clsf == "hard") {
+		    gamma <- t(apply(gamma,1,ind.max))
+		}
 		LL <- -1e10
 		
 		for(i in 1:ns) {

Modified: pkg/depmixS4/R/em.control.R
===================================================================
--- pkg/depmixS4/R/em.control.R	2012-06-12 14:31:19 UTC (rev 517)
+++ pkg/depmixS4/R/em.control.R	2012-06-14 18:46:05 UTC (rev 518)
@@ -1,3 +1,4 @@
 em.control <- 
-function(maxit=500,tol=1e-8,crit="relative",random.start=TRUE) {
-	return(list(maxit=maxit,tol=tol,crit=crit,random.start=random.start))}
\ No newline at end of file
+function(maxit=500,tol=1e-8,crit="relative",random.start=TRUE,classification=c("soft","hard")) {
+    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-12 14:31:19 UTC (rev 517)
+++ pkg/depmixS4/man/em.control.Rd	2012-06-14 18:46:05 UTC (rev 518)
@@ -8,7 +8,7 @@
 
 \usage{
 	
-	em.control(maxit = 500, tol = 1e-08, crit = "relative", random.start = TRUE)
+	em.control(maxit = 500, tol = 1e-08, crit = "relative", random.start = TRUE, classification = c("soft","hard"))
 	
 }
 
@@ -24,6 +24,9 @@
 	\item{random.start}{This is used for a (limited) random
 	initialization of the parameters. See Details.}
 	
+	\item{classification}{Type of classification to states
+	used. See Details.}
+	
 }
 
 \details{
@@ -52,6 +55,14 @@
 given to distinguish between the states.  It is also useful for repeated
 estimation from different starting values.
 
+Argument \code{classification} is used to choose either soft (default) or
+hard classification to states. When using soft classification, observations
+are assigned to states with a weight equal to the posterior probability of
+the state. When using hard classification, observations are assigned to states
+according to the maximum a posteriori (MAP) states (i.e., each observation
+is assigned to one state, which is determined by the Viterbi algorithm in the
+case of \code{depmix} models).
+
 }
 
 \references{
@@ -68,4 +79,4 @@
 
 \author{Ingmar Visser & Maarten Speekenbrink}
 
-\keyword{methods}
\ No newline at end of file
+\keyword{methods}