SCM

SCM Repository

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

View of /trunk/R/fb.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 101 - (download) (annotate)
Sun Mar 23 14:29:34 2008 UTC (11 years, 6 months ago) by ingmarvisser
File size: 2064 byte(s)
Moved latest version to pkg directory
# 
# Maarten Speekenbrink
# 
# FORWARD-BACKWARD algoritme, 23-3-2008
# 

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
}


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