# SCM Repository

[depmix] View of /trunk/R/fb.R
 [depmix] / trunk / R / fb.R

# View of /trunk/R/fb.R

Fri Mar 7 14:40:09 2008 UTC (11 years, 7 months ago) by ingmarvisser
File size: 2027 byte(s)
`Moved generics into allGenerics.R and cleaned up comments etc`
```#
# FORWARD-BACKWARD algoritme
#

fb <- function(init,A,B,ntimes=NULL,return.all=FALSE,stationary=TRUE) {

# Forward-Backward algorithm (used in Baum-Welch)
# Returns alpha, beta, and full data likelihood
# NOTE THE CHANGE IN FROM ROW TO COLUMN SUCH THAT TRANSPOSING A IS NOT NECCESSARY ANYMORE
# IN COMPUTING ALPHA AND BETA BUT IS NOW NECCESSARY IN COMPUTING XI
# A = T*K*K matrix with transition probabilities, from row to column!!!!!!!
# B = T*K matrix with elements ab_{ij} = P(y_i|s_j)
# init = K vector with initial probabilities

# NOTE: to prevent underflow, alpha and beta are scaled, using c

# NOTE: xi[i,j] = P(S[t] = j & S[t+1] = i)

nt <- nrow(B)
ns <- ncol(init)

alpha <- matrix(ncol=ns,nrow=nt)
beta <- matrix(ncol=ns,nrow=nt)
sca <- vector(length=nt)
xi <- array(dim=c(nt,ns,ns))

if(is.null(ntimes)) ntimes <- nt

lt <- length(ntimes)
et <- cumsum(ntimes)
bt <- c(1,et[-lt]+1)

for(case in 1:lt) {
alpha[bt[case],] <- init[case,]*B[bt[case],] # initialize
sca[bt[case]] <- 1/sum(alpha[bt[case],])
alpha[bt[case],] <- alpha[bt[case],]*sca[bt[case]]

if(ntimes[case]>1) {
for(i in bt[case]:(et[case]-1)) {
if(stationary) alpha[i+1,] <- (A[1,,]%*%alpha[i,])*B[i+1,]
else alpha[i+1,] <- (A[i,,]%*%alpha[i,])*B[i+1,]
sca[i+1] <- 1/sum(alpha[i+1,])
alpha[i+1,] <- sca[i+1]*alpha[i+1,]
}
}

beta[et[case],] <- 1*sca[et[case]] # initialize

if(ntimes[case]>1) {
for(i in (et[case]-1):bt[case]) {
if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
else beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[i,,]*sca[i]
}

for(i in bt[case]:(et[case]-1)) {
if(stationary) xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[1,,])
else xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[i,,])
}
}

}

gamma <- alpha*beta/sca
like <- -sum(log(sca))

if(return.all) {
res <- list(alpha=alpha,beta=beta,gamma=gamma,xi=xi,sca=sca,logLike=like)
} else {
res <- list(gamma=gamma,xi=xi,logLike=like)
}

res
}

```