[Depmix-commits] r140 - trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 9 15:08:25 CEST 2008
Author: maarten
Date: 2008-06-09 15:08:25 +0200 (Mon, 09 Jun 2008)
New Revision: 140
Modified:
trunk/R/responseGLM.R
trunk/R/responseGLMBINOM.R
trunk/R/responseGLMGAMMA.R
trunk/R/responseGLMMULTINOM.R
trunk/R/responseGLMPOISSON.R
trunk/R/responseMVN.R
trunk/R/responseNORM.R
trunk/R/simulate.R
Log:
- updated simulate method.
- simulate(response,nsim,...) now returns a (T*nsim)*ncol(response at y) matrix
- simulate(depmix,nsim,...) now returns an object of class depmix.sim with additional slot for states
Modified: trunk/R/responseGLM.R
===================================================================
--- trunk/R/responseGLM.R 2008-06-09 12:01:32 UTC (rev 139)
+++ trunk/R/responseGLM.R 2008-06-09 13:08:25 UTC (rev 140)
@@ -163,7 +163,7 @@
setMethod("fit","GLMresponse",
function(object,w) {
pars <- object at parameters
- fit <- glm.fit(x=object at x,y=object at y,weights=w,family=object at family)
+ fit <- glm.fit(x=object at x,y=object at y,weights=w,family=object at family,start=pars$coefficients)
pars$coefficients <- fit$coefficients
object <- setpars(object,unlist(pars))
object
Modified: trunk/R/responseGLMBINOM.R
===================================================================
--- trunk/R/responseGLMBINOM.R 2008-06-09 12:01:32 UTC (rev 139)
+++ trunk/R/responseGLMBINOM.R 2008-06-09 13:08:25 UTC (rev 140)
@@ -30,7 +30,8 @@
}
nt <- nrow(pr)
response <- rbinom(nt*nsim,size=object at n,prob=pr)
- if(nsim > 1) response <- matrix(response,ncol=nsim)
+ #if(nsim > 1) response <- matrix(response,ncol=nsim)
+ response <- as.matrix(response)
return(response)
}
)
Modified: trunk/R/responseGLMGAMMA.R
===================================================================
--- trunk/R/responseGLMGAMMA.R 2008-06-09 12:01:32 UTC (rev 139)
+++ trunk/R/responseGLMGAMMA.R 2008-06-09 13:08:25 UTC (rev 140)
@@ -29,7 +29,8 @@
}
nt <- nrow(shape)
response <- rgamma(nt*nsim,shape=shape)
- if(nsim > 1) response <- matrix(response,ncol=nsim)
+ # if(nsim > 1) response <- matrix(response,ncol=nsim)
+ response <- as.matrix(response)
return(response)
}
)
Modified: trunk/R/responseGLMMULTINOM.R
===================================================================
--- trunk/R/responseGLMMULTINOM.R 2008-06-09 12:01:32 UTC (rev 139)
+++ trunk/R/responseGLMMULTINOM.R 2008-06-09 13:08:25 UTC (rev 140)
@@ -52,7 +52,8 @@
}
nt <- nrow(pr)
sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nt))
- response <- t(apply(sims,c(2,3), function(x) which(x==1)))
- return(response)
+ sims <- matrix(aperm(sims,c(3,2,1)),nrow=nsim*nt,ncol=ncol(pr))
+ #response <- t(apply(sims,c(2,3), function(x) which(x==1)))
+ return(sims)
}
)
Modified: trunk/R/responseGLMPOISSON.R
===================================================================
--- trunk/R/responseGLMPOISSON.R 2008-06-09 12:01:32 UTC (rev 139)
+++ trunk/R/responseGLMPOISSON.R 2008-06-09 13:08:25 UTC (rev 140)
@@ -29,7 +29,8 @@
}
nt <- nrow(lambda)
response <- rpois(nt*nsim,lambda=lambda)
- if(nsim > 1) response <- matrix(response,ncol=nsim)
+ #if(nsim > 1) response <- matrix(response,ncol=nsim)
+ response <- as.matrix(response)
return(response)
}
)
Modified: trunk/R/responseMVN.R
===================================================================
--- trunk/R/responseMVN.R 2008-06-09 12:01:32 UTC (rev 139)
+++ trunk/R/responseMVN.R 2008-06-09 13:08:25 UTC (rev 140)
@@ -71,7 +71,7 @@
}
nt <- nrow(mu)
response <- mvrnorm(nt*nsim,mu=mu,Sigma=object at parameters$Sigma)
- if(nsim > 1) response <- array(response,dim=c(nt,ncol(response),nsim))
+ #if(nsim > 1) response <- array(response,dim=c(nt,ncol(response),nsim))
return(response)
}
)
Modified: trunk/R/responseNORM.R
===================================================================
--- trunk/R/responseNORM.R 2008-06-09 12:01:32 UTC (rev 139)
+++ trunk/R/responseNORM.R 2008-06-09 13:08:25 UTC (rev 140)
@@ -45,7 +45,8 @@
nt <- length(mu)
sd <- getpars(object)["sd"]
response <- rnorm(nt*nsim,mean=mu,sd=sd)
- if(nsim > 1) response <- matrix(response,ncol=nsim)
+ #if(nsim > 1) response <- matrix(response,ncol=nsim)
+ response <- as.matrix(response)
return(response)
}
)
Modified: trunk/R/simulate.R
===================================================================
--- trunk/R/simulate.R 2008-06-09 12:01:32 UTC (rev 139)
+++ trunk/R/simulate.R 2008-06-09 13:08:25 UTC (rev 140)
@@ -2,6 +2,13 @@
# TODO: move this to all generics etc...
#setGeneric("is.stationary", function(object,...) standardGeneric("is.stationary"))
+setClass("depmix.sim",
+ contains="depmix",
+ representation(
+ states="matrix"
+ )
+)
+
setMethod("is.stationary",signature(object="depmix"),
function(object) {
return(object at stationary)
@@ -35,27 +42,53 @@
states[i,] <- sims[cbind(i,states[i-1,],1:nsim)]
}
}
-# for(j in 1:nsim) {
-# if(is.stationary(object)) {
-# states[i,j] <- simulate(object at transition[[states[i-1,j]]],nsim=1)
-# } else {
-# states[i,j] <- simulate(object at transition[[states[i-1,j]]],nsim=1,time=i)
-# }
-# }
-# }
-# }
- responses <- array(,dim=c(nt,nr,nsim))
+ states <- as.vector(states)
+ responses <- list(length=nr)
+ #responses <- array(,dim=c(nt,nr,nsim))
for(i in 1:nr) {
- tmp <- array(,dim=c(nt,ns,nsim))
+ tmp <- matrix(,nrow=nt*nsim,ncol=NCOL(object at response[[1]][[i]]@y))
for(j in 1:ns) {
- tmp[,j,] <- simulate(object at response[[j]][[i]],nsim=nsim)
+ tmp[states==j,] <- simulate(object at response[[j]][[i]],nsim=nsim)[states==j,]
}
- for(j in 1:nsim) {
- responses[,i,j] <- tmp[cbind(1:nt,states[,j],j)]
+ responses[[i]] <- tmp
+ }
+
+ # generate new depmix.sim object
+ class(object) <- "depmix.sim"
+ object at states <- as.matrix(states)
+
+ object at prior@x <- 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))
}
}
- return(list(states=states,responses=responses))
+ object at ntimes <- rep(object at ntimes,nsim)
+
+ # make appropriate array for transition densities
+ nt <- sum(object at ntimes)
+ if(is.stationary(object)) trDens <- array(0,c(1,ns,ns)) else trDens <- array(0,c(nt,ns,ns))
+
+ # 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
+ }
+ trDens[,,i] <- dens(object at transition[[i]])
+ }
+
+ # compute initial state probabilties
+ object at init <- dens(object at prior)
+ object at trDens <- trDens
+ object at dens <- dns
+
+ return(object)
}
)
More information about the depmix-commits
mailing list