SCM

SCM Repository

[depmix] Annotation of /trunk/R/viterbi.R
ViewVC logotype

Annotation of /trunk/R/viterbi.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 66 - (view) (download)

1 : maarten 55 viterbi <- function(hmm) {
2 : maarten 44 # returns the most likely state sequence
3 : maarten 55 nt <- sum(hmm@ntimes)
4 :     lt <- length(hmm@ntimes)
5 :     et <- cumsum(hmm@ntimes)
6 :     bt <- c(1,et[-lt]+1)
7 :    
8 :     ns <- hmm@nstates
9 :    
10 :     delta <- psi <- matrix(nrow=nt,ncol=ns)
11 :     state <- vector(length=nt)
12 :    
13 :     prior <- exp(logDens(hmm@initModel))
14 :    
15 :     A <- hmm@trans
16 :     B <- apply(hmm@logdens,c(1,3),sum)
17 :    
18 :     for(case in 1:lt) {
19 :     # initialization
20 :     delta[bt[case],] <- - (log(prior[case,]) + B[bt[case],])
21 :     psi[bt[case],] <- 0
22 :     # recursion
23 :     for(i in ((bt[case]+1):et[case])) {
24 :     for(j in 1:ns) {
25 :     delta[i,j] <- min(delta[i-1,] - log(A[i,,j])) - B[i,j]
26 :     k <- which.min(delta[i-1,] - log(A[i,,j]))
27 :     if(length(k) == 0) k <- 0
28 :     psi[i,j] <- k
29 :     }
30 :     }
31 :     #trace maximum likely state
32 :     state[et[case]] <- which.min(delta[et[case],])
33 :     for(i in (et[case]-1):bt[case]) {
34 :     state[i] <- psi[i+1,state[i+1]]
35 :     }
36 :     }
37 :     return(state)
38 :     }
39 :    
40 :     viterbi.fb <- function(A,B,prior) {
41 :     # returns the most likely state sequence
42 : maarten 44 nt <- nrow(B)
43 :     ns <- ncol(A)
44 :     delta <- psi <- matrix(nrow=nt,ncol=ns)
45 :     state <- vector(length=nt)
46 :     # initialization
47 :     delta[1,] <- - (log(prior) + log(B[1,]))
48 :     psi[1,] <- 0
49 :     # recursion
50 :     for(i in 2:nt) {
51 :     for(j in 1:ns) {
52 :     delta[i,j] <- min(delta[i-1,] - log(A[,j])) - log(B[i,j])
53 :     k <- which.min(delta[i-1,] - log(A[,j]))
54 :     if(length(k) == 0) k <- 0
55 :     psi[i,j] <- k
56 :     }
57 :     }
58 :     #trace maximum likely state
59 :     state[nt] <- which.min(delta[nt,])
60 :     for(i in (nt-1):1) {
61 :     state[i] <- psi[i+1,state[i+1]]
62 :     }
63 :     return(state)
64 :     }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge