[Depmix-commits] r224 - trunk/R

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

Author: ingmarvisser
Date: 2008-08-26 15:16:23 +0200 (Tue, 26 Aug 2008)
New Revision: 224

Modified:
   trunk/R/depmix-class.R
   trunk/R/depmixsim-class.R
   trunk/R/responseGLM.R
Log:
Added simulate functie voor mix models

Modified: trunk/R/depmix-class.R
===================================================================
--- trunk/R/depmix-class.R	2008-08-20 14:03:10 UTC (rev 223)
+++ trunk/R/depmix-class.R	2008-08-26 13:16:23 UTC (rev 224)
@@ -44,7 +44,95 @@
 	function(object) return(object at nresp)
 )
 
+setMethod("is.stationary",signature(object="mix"),
+  function(object) {
+		return(TRUE)
+	}
+)
 
+setMethod("simulate",signature(object="mix"),
+	function(object,nsim=1,seed=NULL,...) {
+		
+		if(!is.null(seed)) set.seed(seed)
+		
+		ntim <- ntimes(object)
+		nt <- sum(ntim)
+		bt <- 1:nt
+		
+		nr <- nresp(object)
+		ns <- nstates(object)
+		
+		# simulate state sequences first, then observations
+		
+		# random generation is slow when done separately for each t, so first draw
+		#   variates for all t, and then determine state sequences iteratively
+		states <- array(,dim=c(nt,nsim))
+		states[bt,] <- simulate(object at prior,n=nsim,is.prior=T)
+		sims <- array(,dim=c(nt,ns,nsim))
+		
+# 		for(i in 1:ns) {
+# 			if(is.stationary(object)) {
+# 				# TODO: this is a temporary fix!!! 
+# 				sims[,i,] <- simulate(object at transition[[i]],nsim=nsim,times=rep(1,nt))
+# 			} else {
+# 				sims[,i,] <- simulate(object at transition[[i]],nsim=nsim)
+# 			}
+# 		}
+# 		# track states
+# 		for(case in 1:lt) {
+# 			for(i in (bt[case]+1):et[case]) {
+# 				states[i,] <- sims[cbind(i,states[i-1,],1:nsim)]
+# 			}
+# 		}
+		
+		states <- as.vector(states)
+		responses <- list(length=nr)
+		#responses <- array(,dim=c(nt,nr,nsim))
+		for(i in 1:nr) {
+			tmp <- matrix(,nrow=nt*nsim,ncol=NCOL(object at response[[1]][[i]]@y))
+			for(j in 1:ns) {
+				tmp[states==j,] <- simulate(object at response[[j]][[i]],nsim=nsim)[states==j,]
+			}
+			responses[[i]] <- tmp
+		}
+		
+		# generate new depmix.sim object
+		class(object) <- "mix.sim"
+		object at states <- as.matrix(states)
+		
+		object at prior@x <- as.matrix(apply(object at prior@x,2,rep,nsim))
+		for(j in 1:ns) {
+# 			if(!is.stationary(object)) object at transition[[j]]@x <- as.matrix(apply(object at transition[[j]]@x,2,rep,nsim))
+			for(i in 1:nr) {
+				object at response[[j]][[i]]@y <- as.matrix(responses[[i]])
+				object at response[[j]][[i]]@x <- as.matrix(apply(object at response[[j]][[i]]@x,2,rep,nsim))
+			}
+		}
+		object at ntimes <- rep(object at ntimes,nsim)
+		
+		# make appropriate array for transition densities
+		nt <- sum(object at ntimes)
+		
+		# make appropriate array for response densities
+		dns <- array(,c(nt,nr,ns))
+		
+		# compute observation and transition densities
+		for(i in 1:ns) {
+			for(j in 1:nr) {
+				dns[,j,i] <- dens(object at response[[i]][[j]]) # remove this response as an argument from the call to setpars
+			}
+		}
+		
+		# compute initial state probabilties
+		object at init <- dens(object at prior)
+		object at dens <- dns
+		
+		return(object)
+	}
+)
+
+
+
 # 
 # PRINT method
 # 

Modified: trunk/R/depmixsim-class.R
===================================================================
--- trunk/R/depmixsim-class.R	2008-08-20 14:03:10 UTC (rev 223)
+++ trunk/R/depmixsim-class.R	2008-08-26 13:16:23 UTC (rev 224)
@@ -5,4 +5,11 @@
   representation(
     states="matrix"
   )
+)
+
+setClass("mix.sim",
+  contains="mix",
+  representation(
+	states="matrix"
+  )
 )
\ No newline at end of file

Modified: trunk/R/responseGLM.R
===================================================================
--- trunk/R/responseGLM.R	2008-08-20 14:03:10 UTC (rev 223)
+++ trunk/R/responseGLM.R	2008-08-26 13:16:23 UTC (rev 224)
@@ -60,7 +60,8 @@
 		}
 		if(family$family=="multinomial") {
 			y <- model.response(mf)
-			if(is.factor(y)) y <- model.matrix(~y-1) else if(is.numeric(y)) y <- model.matrix(~factor(y)-1)
+			# FIX THIS: comment out when it produces errors when the response is a matrix, see multinom help page
+ 			if(is.factor(y)) y <- model.matrix(~y-1) else if(is.numeric(y)) y <- model.matrix(~factor(y)-1)
 			parameters$coefficients <- matrix(0,ncol=ncol(y),nrow=ncol(x))
 			if(is.null(fixed)) {
 				fixed <- parameters$coefficients