[Depmix-commits] r216 - trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 1 13:56:48 CEST 2008
Author: ingmarvisser
Date: 2008-08-01 13:56:47 +0200 (Fri, 01 Aug 2008)
New Revision: 216
Modified:
trunk/R/depmix-class.R
trunk/R/depmixfit.R
trunk/R/viterbi.R
Log:
Fixed bug in viterbi.
Modified: trunk/R/depmix-class.R
===================================================================
--- trunk/R/depmix-class.R 2008-07-29 14:35:28 UTC (rev 215)
+++ trunk/R/depmix-class.R 2008-08-01 11:56:47 UTC (rev 216)
@@ -125,86 +125,88 @@
)
setMethod("simulate",signature(object="depmix"),
- function(object,nsim=1,seed=NULL,...) {
- if(!is.null(seed)) set.seed(seed)
- 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) {
- 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) <- "depmix.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)
- 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,...) {
+
+ if(!is.null(seed)) set.seed(seed)
+
+ 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) {
+ 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) <- "depmix.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)
+ 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: trunk/R/depmixfit.R
===================================================================
--- trunk/R/depmixfit.R 2008-07-29 14:35:28 UTC (rev 215)
+++ trunk/R/depmixfit.R 2008-08-01 11:56:47 UTC (rev 216)
@@ -84,15 +84,11 @@
lin.l <- c(lin.l,conrows.lower)
}
}
-
- print(lincon)
-
+
# select only those columns of the constraint matrix that correspond to non-fixed parameters
linconFull <- lincon
lincon <- lincon[,!fixed,drop=FALSE]
-
- print(lincon)
-
+
# set donlp2 control parameters
cntrl <- donlp2.control(hessian=FALSE,difftype=2,report=TRUE)
Modified: trunk/R/viterbi.R
===================================================================
--- trunk/R/viterbi.R 2008-07-29 14:35:28 UTC (rev 215)
+++ trunk/R/viterbi.R 2008-08-01 11:56:47 UTC (rev 216)
@@ -27,26 +27,27 @@
psi[bt[case],] <- 0
# recursion
if(object at ntimes[case]>1) {
- for(i in ((bt[case]+1):et[case])) {
+ for(tt in ((bt[case]+1):et[case])) {
for(j in 1:ns) {
if(!object at stationary) {
- delta[i,j] <- max(delta[i-1,]*(A[i,,j]))*B[i,j]
- k <- which.max(delta[i-1,]*A[i,,j])
+ delta[tt,j] <- max(delta[tt-1,]*(A[tt,j,]))*B[tt,j]
+ k <- which.max(delta[tt-1,]*A[tt,j,])
} else {
- delta[i,j] <- max(delta[i-1,]*A[,,j])*B[i,j]
- k <- which.max(delta[i-1,]*A[,,j])
+ delta[tt,j] <- max(delta[tt-1,]*(A[1,j,]))*B[tt,j]
+ k <- which.max(delta[tt-1,]*A[1,j,])
}
- if(length(k) == 0) k <- 0
- psi[i,j] <- k
+ if(length(k) == 0) k <- 0 # what's this doing here??? can this ever occur? FIX ME
+ psi[tt,j] <- k
}
- delta[i,] <- delta[i,]/(sum(delta[i,]))
+ delta[tt,] <- delta[tt,]/(sum(delta[tt,]))
+
}
}
# trace maximum likely state
state[et[case]] <- which.max(delta[et[case],])
- # this doesn't need a for loop does it???? FIX ME
+ # this doesn't need a for loop does it???? FIX ME
if(object at ntimes[case]>1) {
for(i in (et[case]-1):bt[case]) {
state[i] <- psi[i+1,state[i+1]]
More information about the depmix-commits
mailing list