[Depmix-commits] r439 - trunk

Mon Aug 9 14:43:58 CEST 2010

Author: ingmarvisser
Date: 2010-08-09 14:43:58 +0200 (Mon, 09 Aug 2010)
New Revision: 439

Modified:
   trunk/depmixNew-test7.R
   trunk/start.R
Log:


Modified: trunk/depmixNew-test7.R
===================================================================

--- trunk/depmixNew-test7.R	2010-08-09 12:43:15 UTC (rev 438)
+++ trunk/depmixNew-test7.R	2010-08-09 12:43:58 UTC (rev 439)
@@ -44,28 +44,67 @@
 
 fm2 <- fit(mod,meth="rsolnp")
 
-# Bovenstaande geeft warnings en uiteindelijk error in qr.default() ...
 
 
 
 
-fm2 <- fit(mod,meth="donlp")
+data(speed)
 
+set.seed(1)
+tr=runif(6)
+trst=c(tr[1:2],0,tr[3:5],0,tr[6])
+mod1 <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,nstates=2,
+	family=list(gaussian(),multinomial("identity")),trstart=trst)
+# fit the model
+fmod1 <- fit(mod1)
+fmod1 # to see the logLik and optimization information
+# to see the parameters
+summary(fmod1)
 
+# FIX SOME PARAMETERS
 
+pars <- c(unlist(getpars(fmod1)))
 
+pars[1]=0
+pars[2]=1 # this means the process will always start in state 2
+pars[13]=0.5
+pars[14]=0.5 # the corr parameters in state 1 are now both 0, corresponding the 0.5 prob
+mod2 <- setpars(mod1,pars)
 
+# fix the parameters by setting: 
+free <- c(0,0,rep(c(0,1),4),1,1,0,0,1,1,1,1)
+# fit the model
+fmod2 <- fit(mod2,fixed=!free)
 
+# likelihood ratio insignificant, hence fmod2 better than fmod1
+llratio(fmod1,fmod2)
 
+# NOW ADD SOME GENERAL LINEAR CONSTRAINTs
 
+pars <- c(unlist(getpars(fmod2)))
+pars[4] <- pars[8] <- -4
+pars[6] <- pars[10] <- 10
+mod3 <- setpars(mod2,pars)
 
+logLik(mod3)
 
+# start with fixed and free parameters
+conpat <- c(0,0,rep(c(0,1),4),1,1,0,0,1,1,1,1)
+# constrain the beta's on the transition parameters to be equal
+conpat[4] <- conpat[8] <- 2
+conpat[6] <- conpat[10] <- 3
 
+fmod3 <- fit(mod3,equal=conpat)
 
 
 
 
 
+
+
+
+
+
 rsolnp control list:
 
 rho

Modified: trunk/start.R
===================================================================
--- trunk/start.R	2010-08-09 12:43:15 UTC (rev 438)
+++ trunk/start.R	2010-08-09 12:43:58 UTC (rev 439)
@@ -1,3 +1,131 @@
+
+
+library(depmixS4)
+
+data(balance)
+# four binary items on the balance scale task
+
+instart=c(0.5,0.5)
+set.seed(1)
+respstart=runif(16)
+# note that ntimes argument is used to make this a mixture model
+mod4 <- mix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
+	family=list(multinomial("identity"),multinomial("identity"),
+		multinomial("identity"),multinomial("identity")),
+	respstart=respstart,instart=instart)
+
+fmod4 <- fit(mod4, em.control=list(random.start=TRUE))
+
+fmod4
+
+mod4b <- mix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
+	family=list(multinomial("identity"),multinomial("identity"),
+		multinomial("identity"),multinomial("identity")))
+
+
+fmod4b <- fit(mod4b)
+
+fmod4b <- fit(mod4b, em.control=list(tol=1e-8,
+		crit=c("relative","absolute"),random.start=TRUE))
+
+
+
+# add age as covariate on class membership by using the prior argument
+instart=c(0.5,0.5,0,0) # we need the initial probs and the coefficients of age 
+set.seed(2)
+respstart=runif(16)
+mod5 <- mix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
+	family=list(multinomial("identity"),multinomial("identity"),
+		multinomial("identity"),multinomial("identity")),
+	instart=instart, respstart=respstart, prior=~age, initdata=balance)
+
+fmod5 <- fit(mod5)
+
+mod5b <- mix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
+	family=list(multinomial("identity"),multinomial("identity"),
+		multinomial("identity"),multinomial("identity")),
+	prior=~age, initdata=balance)
+
+fmod5b <- fit(mod5b)
+
+fmod5b <- fit(mod5b,em.control=list(tol=1e-8,
+		crit=c("relative","absolute"),random.start=TRUE))
+
+
+data(speed)
+
+# 2-state model on rt and corr from speed data set 
+# with Pacc as covariate on the transition matrix
+# starting values for the transition pars (without 
+# those EM does not get off the ground)
+set.seed(1)
+tr=runif(6)
+trst=c(tr[1:2],0,tr[3:5],0,tr[6])
+mod1 <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,nstates=2,
+	family=list(gaussian(),multinomial("identity")),trstart=trst)
+# fit the model
+fmod1 <- fit(mod1)
+fmod1 # to see the logLik and optimization information
+# to see the parameters
+
+mod1b <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,nstates=2,
+	family=list(gaussian(),multinomial("identity")))
+# fit the model
+fmod1 <- fit(mod1b)
+
+fmod1 <- fit(mod1b,em.control=list(tol=1e-8,
+		crit=c("relative","absolute"),random.start=TRUE))
+
+
+# generate data from two different multivariate normal distributions
+m1 <- c(0,1)
+sd1 <- matrix(c(1,0.7,.7,1),2,2)
+m2 <- c(1,0)
+sd2 <- matrix(c(2,.1,.1,1),2,2)
+y1 <- mvrnorm(50,m1,sd1)
+y2 <- mvrnorm(50,m2,sd2)
+y <- rbind(y1,y2)
+
+# now use makeDepmix to create a depmix model for this bivariate normal timeseries
+rModels <-  list()
+rModels[[1]] <- list(MVNresponse(y~1))
+rModels[[2]] <- list(MVNresponse(y~1))
+
+trstart=c(0.9,0.1,0.1,0.9)
+
+transition <- list()
+transition[[1]] <- transInit(~1,nstates=2,data=data.frame(1),pstart=c(trstart[1:2]))
+transition[[2]] <- transInit(~1,nstates=2,data=data.frame(1),pstart=c(trstart[3:4]))
+
+instart=runif(2)
+inMod <- transInit(~1,ns=2,ps=instart,data=data.frame(1))
+
+mod <- makeDepmix(response=rModels,transition=transition,prior=inMod)
+
+fm2 <- fit(mod)
+
+
+# now use makeDepmix to create a depmix model for this bivariate normal timeseries
+rModels <-  list()
+rModels[[1]] <- list(MVNresponse(y~1))
+rModels[[2]] <- list(MVNresponse(y~1))
+
+transition <- list()
+transition[[1]] <- transInit(~1,nstates=2,data=data.frame(1))
+transition[[2]] <- transInit(~1,nstates=2,data=data.frame(1))
+
+inMod <- transInit(~1,ns=2,data=data.frame(1))
+
+mod2b <- makeDepmix(response=rModels,transition=transition,prior=inMod)
+
+fm2b <- fit(mod2b,em.control=list(tol=1e-8,
+		crit=c("relative","absolute"),random.start=TRUE))
+
+
+
+
+
+
 # 
 # Maarten Speekenbrink 23-3-2008
 # 

