[Depmix-commits] r97 - pkg trunk trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 21 17:13:42 CET 2008
Author: ingmarvisser
Date: 2008-03-21 17:13:42 +0100 (Fri, 21 Mar 2008)
New Revision: 97
Added:
trunk/depmixNew-test4.R
Modified:
pkg/NAMESPACE
trunk/R/EM.R
trunk/R/depmix.fitted.R
trunk/R/viterbi.R
trunk/depmixNew-test2.R
Log:
Added example file for viterbi and added viterbi2 with different scaling
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2008-03-20 22:30:09 UTC (rev 96)
+++ pkg/NAMESPACE 2008-03-21 16:13:42 UTC (rev 97)
@@ -1,3 +1,4 @@
+
import(methods)
importFrom(stats, predict)
@@ -28,6 +29,7 @@
BIC,
fit,
npar,
+ freepars,
getdf,
nobs,
depmix,
Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R 2008-03-20 22:30:09 UTC (rev 96)
+++ trunk/R/EM.R 2008-03-21 16:13:42 UTC (rev 97)
@@ -79,7 +79,7 @@
object at conMat <- matrix()
# what do we want in slot posterior?
- object at posterior <- fbo$gamma
+ object at posterior <- viterbi(object)
object
}
Modified: trunk/R/depmix.fitted.R
===================================================================
--- trunk/R/depmix.fitted.R 2008-03-20 22:30:09 UTC (rev 96)
+++ trunk/R/depmix.fitted.R 2008-03-21 16:13:42 UTC (rev 97)
@@ -10,7 +10,7 @@
setClass("depmix.fitted",
representation(message="character", # convergence information
conMat="matrix", # constraint matrix on the parameters for general linear constraints
- posterior="matrix" # posterior probabilities for the states
+ posterior="data.frame" # posterior probabilities for the states
),
contains="depmix"
)
@@ -18,7 +18,7 @@
setMethod("fit",
signature(object="depmix"),
function(object,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,...) {
-
+
# when there are linear constraints donlp should be used
# otherwise EM is good
@@ -131,12 +131,10 @@
# put the result back into the model
allpars[!fixed] <- result$par
object <- setpars(object,allpars)
+ object at posterior <- viterbi(object)
- }
+ }
- # FIX ME
- object at posterior <- matrix(0,1,1)
-
return(object)
}
)
Modified: trunk/R/viterbi.R
===================================================================
--- trunk/R/viterbi.R 2008-03-20 22:30:09 UTC (rev 96)
+++ trunk/R/viterbi.R 2008-03-21 16:13:42 UTC (rev 97)
@@ -1,47 +1,112 @@
viterbi <- function(object) {
- # returns the most likely state sequence
- nt <- sum(object at ntimes)
- lt <- length(object at ntimes)
- et <- cumsum(object at ntimes)
- bt <- c(1,et[-lt]+1)
+ # returns the most likely state sequence
+ nt <- sum(object at ntimes)
+ lt <- length(object at ntimes)
+ et <- cumsum(object at ntimes)
+ bt <- c(1,et[-lt]+1)
- ns <- object at nstates
-
- delta <- psi <- matrix(nrow=nt,ncol=ns)
- state <- vector(length=nt)
-
- prior <- object at init
-
+ ns <- object at nstates
+
+ delta <- psi <- matrix(nrow=nt,ncol=ns)
+ state <- vector(length=nt)
+
+ prior <- object at init
+
A <- object at trDens
- B <- apply(log(object at dens),c(1,3),sum)
-
- for(case in 1:lt) {
- # initialization
- delta[bt[case],] <- - (log(prior[case,]) + B[bt[case],])
- psi[bt[case],] <- 0
- # recursion
- for(i in ((bt[case]+1):et[case])) {
- for(j in 1:ns) {
- if(!object at stationary) {
- delta[i,j] <- min(delta[i-1,] - log(A[i,,j])) - B[i,j]
- k <- which.min(delta[i-1,] - log(A[i,,j]))
- } else {
- delta[i,j] <- min(delta[i-1,] - log(A[,,j])) - B[i,j]
- k <- which.min(delta[i-1,] - log(A[,,j]))
- }
- if(length(k) == 0) k <- 0
- psi[i,j] <- k
- }
- }
- #trace maximum likely state
- state[et[case]] <- which.min(delta[et[case],])
- for(i in (et[case]-1):bt[case]) {
- state[i] <- psi[i+1,state[i+1]]
- }
- }
- return(state)
+ B <- apply(log(object at dens),c(1,3),sum)
+
+ for(case in 1:lt) {
+ # initialization
+ delta[bt[case],] <- - (log(prior[case,]) + B[bt[case],])
+ psi[bt[case],] <- 0
+ # recursion
+ if(object at ntimes[case]>1) {
+ for(i in ((bt[case]+1):et[case])) {
+ for(j in 1:ns) {
+ if(!object at stationary) {
+ delta[i,j] <- min(delta[i-1,] - log(A[i,,j])) - B[i,j]
+ k <- which.min(delta[i-1,] - log(A[i,,j]))
+ } else {
+ delta[i,j] <- min(delta[i-1,] - log(A[,,j])) - B[i,j]
+ k <- which.min(delta[i-1,] - log(A[,,j]))
+ }
+ if(length(k) == 0) k <- 0
+ psi[i,j] <- k
+ }
+ }
+ }
+
+ # trace maximum likely state
+ state[et[case]] <- which.min(delta[et[case],])
+
+ # 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]]
+ }
+ }
+ }
+
+ delta <- data.frame(state,delta)
+ return(delta)
}
+viterbi2 <- function(object) {
+ # returns the most likely state sequence
+ nt <- sum(object at ntimes)
+ lt <- length(object at ntimes)
+ et <- cumsum(object at ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ ns <- object at nstates
+
+ delta <- psi <- matrix(nrow=nt,ncol=ns)
+ state <- vector(length=nt)
+
+ prior <- object at init
+
+ A <- object at trDens
+ B <- apply((object at dens),c(1,3),prod)
+
+ for(case in 1:lt) {
+ # initialization
+ delta[bt[case],] <- prior[case,]*B[bt[case],]
+ delta[bt[case],] <- delta[bt[case],]/(sum(delta[bt[case],]))
+ psi[bt[case],] <- 0
+ # recursion
+ if(object at ntimes[case]>1) {
+ for(i 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])
+ } else {
+ delta[i,j] <- max(delta[i-1,]*A[,,j])*B[i,j]
+ k <- which.max(delta[i-1,]*A[,,j])
+ }
+ if(length(k) == 0) k <- 0
+ psi[i,j] <- k
+ }
+ delta[i,] <- delta[i,]/(sum(delta[i,]))
+ }
+ }
+
+ # trace maximum likely state
+ state[et[case]] <- which.max(delta[et[case],])
+
+ # 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]]
+ }
+ }
+ }
+
+ delta <- data.frame(state,delta)
+ return(delta)
+}
+
+
viterbi.fb <- function(A,B,prior) {
# returns the most likely state sequence
nt <- nrow(B)
Modified: trunk/depmixNew-test2.R
===================================================================
--- trunk/depmixNew-test2.R 2008-03-20 22:30:09 UTC (rev 96)
+++ trunk/depmixNew-test2.R 2008-03-21 16:13:42 UTC (rev 97)
@@ -10,6 +10,8 @@
# test model with EM optimization, no covariates
#
+data(speed)
+
trstart=c(0.899,0.101,0.084,0.916)
instart=c(0.5,0.5)
resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
@@ -22,7 +24,17 @@
mod1 <- fit(mod)
ll <- logLik(mod1)
+source("R/depmix.fitted.R")
+source("R/viterbi.R")
+x1=viterbi(mod1)
+x2=viterbi2(mod1)
+round(head(x1),3)
+round(head(x2),3)
+cor(x1[,1],x2[,1])
+
+
+
#
# Test optimization using Rdonlp2
#
Added: trunk/depmixNew-test4.R
===================================================================
--- trunk/depmixNew-test4.R (rev 0)
+++ trunk/depmixNew-test4.R 2008-03-21 16:13:42 UTC (rev 97)
@@ -0,0 +1,37 @@
+
+#
+# Started by Ingmar Visser 26-2-2008
+#
+# Usage: go to trunk directory and source this file in R
+#
+
+#
+# BALANCE SCALE data example with age as covariate on class membership
+#
+
+library(depmixS4)
+
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
+
+
+data(balance)
+
+# 'log Lik.' -1083.036 (df=9)
+
+#
+# Add age as covariate on class membership
+#
+
+instart=c(0.5,0.5,0,0)
+respstart=c(rep(c(0.1,0.9),4),rep(c(0.9,0.1),4))
+trstart=c(1,0,0,1)
+mod2 <- depmix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
+ family=list(multinomial(),multinomial(),multinomial(),multinomial()),
+ trstart=trstart, instart=instart, respstart=respstart,
+ ntimes=rep(1,nrow(balance)), prior=~age, initdata=balance)
+
+fixed <- c(1,0,1,0,1,1,1,1,rep(c(1,0),8))
+mod3 <- fit(mod2,fixed=fixed)
+
+
+
Property changes on: trunk/depmixNew-test4.R
___________________________________________________________________
Name: svn:executable
+ *
More information about the depmix-commits
mailing list