[Depmix-commits] r212 - pkg/R pkg/man trunk/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 2 00:07:39 CEST 2008
Author: ingmarvisser
Date: 2008-07-02 00:07:39 +0200 (Wed, 02 Jul 2008)
New Revision: 212
Modified:
pkg/R/depmix-class.R
pkg/man/simulate.Rd
trunk/man/simulate.Rd
Log:
Added example to simulate (which does not work now)
Modified: pkg/R/depmix-class.R
===================================================================
--- pkg/R/depmix-class.R 2008-07-01 21:01:20 UTC (rev 211)
+++ pkg/R/depmix-class.R 2008-07-01 22:07:39 UTC (rev 212)
@@ -125,80 +125,81 @@
)
setMethod("simulate",signature(object="depmix"),
- function(object,nsim=1,seed=NULL,...) {
- ntim <- ntimes(object)
- nt <- sum(ntim)
- lt <- length(ntim)
- et <- cumsum(ntim)
- bt <- c(1,et[-lt]+1)
-
- 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) {
- 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) <- "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))
- }
- }
- 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)
- }
+ function(object,nsim=1,seed=NULL,...) {
+
+ ntim <- ntimes(object)
+ nt <- sum(ntim)
+ lt <- length(ntim)
+ et <- cumsum(ntim)
+ bt <- c(1,et[-lt]+1)
+
+ 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) {
+ 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) <- "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))
+ }
+ }
+ 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)
+ }
)
Modified: pkg/man/simulate.Rd
===================================================================
--- pkg/man/simulate.Rd 2008-07-01 21:01:20 UTC (rev 211)
+++ pkg/man/simulate.Rd 2008-07-01 22:07:39 UTC (rev 212)
@@ -78,6 +78,20 @@
}
+\examples{
+
+y <- rnorm(1000)
+respst <- c(0,1,2,1)
+trst <- c(0.9,0.1,0.1,0.9)
+
+df <- data.frame(y=y)
+
+mod <- depmix(y~1,data=df,respst=respst,trst=trst,inst=c(0.5,0.5),nti=1000,nst=2)
+
+mod <- simulate(mod)
+
+}
+
\author{Maarten Speekenbrink}
\keyword{methods}
Modified: trunk/man/simulate.Rd
===================================================================
--- trunk/man/simulate.Rd 2008-07-01 21:01:20 UTC (rev 211)
+++ trunk/man/simulate.Rd 2008-07-01 22:07:39 UTC (rev 212)
@@ -78,6 +78,19 @@
}
+\examples{
+
+y <- rnorm(1000)
+respst <- c(0,1,2,1)
+trst <- c(0.9,0.1,0.1,0.9)
+
+df <- data.frame(y=y)
+
+mod <- depmix(y~1,data=df,respst=respst,trst=trst,inst=c(0.5,0.5),nti=1000,nst=2)
+
+mod <- simulate(mod)
+
+}
\author{Maarten Speekenbrink}
\keyword{methods}
More information about the depmix-commits
mailing list