[Depmix-commits] r223 - in trunk: . tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 20 16:03:10 CEST 2008
Author: ingmarvisser
Date: 2008-08-20 16:03:10 +0200 (Wed, 20 Aug 2008)
New Revision: 223
Added:
trunk/tests/test1speed.R
trunk/tests/test1speed.Rout.save
trunk/tests/test2getsetpars.R
trunk/tests/test2getsetpars.Rout.save
trunk/tests/test3responses.R
Removed:
trunk/tests/depmixNew-test1.R
Modified:
trunk/NEWS
Log:
Restructured test files; added tests for response models; added .Rout.save files so errors are checked in R CMD check
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2008-08-20 13:14:18 UTC (rev 222)
+++ trunk/NEWS 2008-08-20 14:03:10 UTC (rev 223)
@@ -3,8 +3,11 @@
o Fixed a bug in the Viterbi algorithm used to compute posterior states
+ o Restructured test files somewhat
+
o
+
Changes in depmixS4 version 0.2-0
o restructured R and Rd (help) files; added depmixS4 help with a short
Deleted: trunk/tests/depmixNew-test1.R
===================================================================
--- trunk/tests/depmixNew-test1.R 2008-08-20 13:14:18 UTC (rev 222)
+++ trunk/tests/depmixNew-test1.R 2008-08-20 14:03:10 UTC (rev 223)
@@ -1,151 +0,0 @@
-
-#
-# Started by Ingmar Visser & Maarten Speekenbrink, 31-1-2007
-#
-# Usage: go to trunk directory and source this file in R, if the program
-# still works it should return TRUE at every test (or make immediate sense
-# otherwise)
-
-# Changes:
-#
-# 16-5-2007, made this a new copy with basic tests of computing the
-# likelihood, copying parameters etc
-#
-# Other tests with optimization of models are moved to depmix-test2.R
-#
-
-#
-# TEST 1: speed data model with optimal parameters, compute the likelihood
-#
-
-require(depmixS4)
-
-data(speed)
-
-pars <- c(1,0.916,0.084,0.101,0.899,6.39,0.24,0.098,0.902,5.52,0.202,0.472,0.528,1,0)
-
-rModels <- list(
- list(
- GLMresponse(formula=rt~1,data=speed,family=gaussian(),pstart=c(5.52,.202)),
- GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(0.472,0.528))
- ),
- list(
- GLMresponse(formula=rt~1,data=speed,family=gaussian(),pstart=c(6.39,.24)),
- GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(.098,.902))
- )
-)
-
-trstart=c(0.899,0.101,0.084,0.916)
-
-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=c(0,1)
-inMod <- transInit(~1,ns=2,ps=instart,data=data.frame(rep(1,3)))
-
-mod <- makeDepmix(response=rModels,transition=transition,prior=inMod,ntimes=attr(speed,"ntimes"))
-
-ll <- logLik(mod)
-ll.fb <- logLik(mod,method="fb")
-
-logl <- -296.115107102 # see above
-
-cat("Test 1: ", all.equal(c(ll),logl,check.att=FALSE), "(loglike of speed data) \n")
-
-
-#
-# model specification made easy
-#
-
-library(depmixS4)
-
-resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
-trstart=c(0.899,0.101,0.084,0.916)
-instart=c(0,1)
-mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()),respstart=resp,trstart=trstart,instart=instart,prob=T)
-
-ll2 <- logLik(mod)
-
-cat("Test 1b: ", all.equal(c(ll),c(ll2),check.att=FALSE), "(loglike of speed data) \n")
-
-#
-# TEST 2
-#
-# To check the density function for the multinomial responses with a covariate
-# test a model with a single state, which should be identical to a glm
-# first fit a model without covariate
-#
-
-invlogit <- function(lp) {exp(lp)/(1+exp(lp))}
-
-acc <- glm(corr~1,data=speed,family=binomial)
-
-p1 <- invlogit(coef(acc)[1])
-p0 <- 1-p1
-
-mod <- depmix(corr~1,data=speed,nst=1,family=multinomial(),trstart=1,instart=c(1),respstart=c(p0,p1))
-
-ll <- logLik(mod)
-
-dev <- -2*ll
-
-cat("Test 2: ", all.equal(c(dev),acc$deviance),"(loglike of 1-comp glm on acc data) \n")
-
-
-#
-# TEST 3
-#
-# now add the covariate and compute the loglikelihood
-#
-
-acc <- glm(corr~Pacc,data=speed,family=binomial)
-
-p1 <- invlogit(coef(acc)[1])
-p0 <- 1-p1
-
-pstart=c(p0,p1,0,coef(acc)[2])
-
-mod <- depmix(corr~Pacc,data=speed,family=multinomial(),trstart=1,instart=1,respst=pstart,nstate=1)
-
-ll <- logLik(mod)
-dev <- -2*ll
-
-cat("Test 3: ", all.equal(c(dev),acc$deviance),"(same but now with covariate) \n")
-
-#
-# 2-state model with covariate
-#
-
-trstart=c(0.896,0.104,0.084,0.916)
-trstart=c(trstart[1:2],0,0.01,trstart[3:4],0,0.01)
-instart=c(0,1)
-resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
-
-mod <- depmix(list(rt~1,corr~1),data=speed,family=list(gaussian(),multinomial()),transition=~Pacc,trstart=trstart,instart=instart,respst=resp,nst=2)
-ll <- logLik(mod)
-
-cat("Test 4: ll is now larger than speedll, ie ll is better due to introduction of a covariate \n")
-cat("Test 4: ", ll,"\t", logl, "\n")
-cat("Test 4: ", ll > logl, "\n")
-
-
-#
-# test getpars and setpars
-#
-
-trstart=c(0.896,0.104,0.084,0.916)
-trstart=c(trstart[1:2],0,0.01,trstart[3:4],0,0.01)
-instart=c(0,1)
-resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
-
-mod <- depmix(list(rt~1,corr~1),data=speed,family=list(gaussian(),multinomial()),transition=~Pacc,trstart=trstart,instart=instart,respst=resp,nst=2)
-
-mod1 <- setpars(mod,getpars(mod))
-
-cat("Test 5: ", all.equal(getpars(mod),getpars(mod1)), "(getpars and setpars) \n")
-
-
-
-
-
Added: trunk/tests/test1speed.R
===================================================================
--- trunk/tests/test1speed.R (rev 0)
+++ trunk/tests/test1speed.R 2008-08-20 14:03:10 UTC (rev 223)
@@ -0,0 +1,114 @@
+#
+# TEST 1: speed data model with optimal parameters, compute the likelihood
+#
+
+require(depmixS4)
+
+data(speed)
+
+pars <- c(1,0.916,0.084,0.101,0.899,6.39,0.24,0.098,0.902,5.52,0.202,0.472,0.528,1,0)
+
+rModels <- list(
+ list(
+ GLMresponse(formula=rt~1,data=speed,family=gaussian(),pstart=c(5.52,.202)),
+ GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(0.472,0.528))
+ ),
+ list(
+ GLMresponse(formula=rt~1,data=speed,family=gaussian(),pstart=c(6.39,.24)),
+ GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(.098,.902))
+ )
+)
+
+trstart=c(0.899,0.101,0.084,0.916)
+
+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=c(0,1)
+inMod <- transInit(~1,ns=2,ps=instart,data=data.frame(rep(1,3)))
+
+mod <- makeDepmix(response=rModels,transition=transition,prior=inMod,ntimes=attr(speed,"ntimes"))
+
+ll <- logLik(mod)
+ll.fb <- logLik(mod,method="fb")
+
+logl <- -296.115107102 # see above
+
+cat("Test 1: ", all.equal(c(ll),logl,check.att=FALSE), "(loglike of speed data) \n")
+
+
+#
+# model specification made easy
+#
+
+library(depmixS4)
+
+resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+trstart=c(0.899,0.101,0.084,0.916)
+instart=c(0,1)
+mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()),respstart=resp,trstart=trstart,instart=instart,prob=T)
+
+ll2 <- logLik(mod)
+
+cat("Test 1b: ", all.equal(c(ll),c(ll2),check.att=FALSE), "(loglike of speed data) \n")
+
+#
+# TEST 2
+#
+# To check the density function for the multinomial responses with a covariate
+# test a model with a single state, which should be identical to a glm
+# first fit a model without covariate
+#
+
+invlogit <- function(lp) {exp(lp)/(1+exp(lp))}
+
+acc <- glm(corr~1,data=speed,family=binomial)
+
+p1 <- invlogit(coef(acc)[1])
+p0 <- 1-p1
+
+mod <- depmix(corr~1,data=speed,nst=1,family=multinomial(),trstart=1,instart=c(1),respstart=c(p0,p1))
+
+ll <- logLik(mod)
+
+dev <- -2*ll
+
+cat("Test 2: ", all.equal(c(dev),acc$deviance),"(loglike of 1-comp glm on acc data) \n")
+
+
+#
+# TEST 3
+#
+# now add the covariate and compute the loglikelihood
+#
+
+acc <- glm(corr~Pacc,data=speed,family=binomial)
+
+p1 <- invlogit(coef(acc)[1])
+p0 <- 1-p1
+
+pstart=c(p0,p1,0,coef(acc)[2])
+
+mod <- depmix(corr~Pacc,data=speed,family=multinomial(),trstart=1,instart=1,respst=pstart,nstate=1)
+
+ll <- logLik(mod)
+dev <- -2*ll
+
+cat("Test 3: ", all.equal(c(dev),acc$deviance),"(same but now with covariate) \n")
+
+#
+# 2-state model with covariate
+#
+
+trstart=c(0.896,0.104,0.084,0.916)
+trstart=c(trstart[1:2],0,0.01,trstart[3:4],0,0.01)
+instart=c(0,1)
+resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+
+mod <- depmix(list(rt~1,corr~1),data=speed,family=list(gaussian(),multinomial()),transition=~Pacc,trstart=trstart,instart=instart,respst=resp,nst=2)
+ll <- logLik(mod)
+
+cat("Test 4: ll is now larger than speedll, ie ll is better due to introduction of a covariate \n")
+cat("Test 4: ", ll,"\t", logl, "\n")
+cat("Test 4: ", ll > logl, "\n")
Property changes on: trunk/tests/test1speed.R
___________________________________________________________________
Name: svn:executable
+ *
Added: trunk/tests/test1speed.Rout.save
===================================================================
--- trunk/tests/test1speed.Rout.save (rev 0)
+++ trunk/tests/test1speed.Rout.save 2008-08-20 14:03:10 UTC (rev 223)
@@ -0,0 +1,141 @@
+
+R version 2.7.1 (2008-06-23)
+Copyright (C) 2008 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+ Natural language support but running in an English locale
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> #
+> # TEST 1: speed data model with optimal parameters, compute the likelihood
+> #
+>
+> require(depmixS4)
+>
+> data(speed)
+>
+> pars <- c(1,0.916,0.084,0.101,0.899,6.39,0.24,0.098,0.902,5.52,0.202,0.472,0.528,1,0)
+>
+> rModels <- list(
++ list(
++ GLMresponse(formula=rt~1,data=speed,family=gaussian(),pstart=c(5.52,.202)),
++ GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(0.472,0.528))
++ ),
++ list(
++ GLMresponse(formula=rt~1,data=speed,family=gaussian(),pstart=c(6.39,.24)),
++ GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(.098,.902))
++ )
++ )
+>
+> trstart=c(0.899,0.101,0.084,0.916)
+>
+> 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=c(0,1)
+> inMod <- transInit(~1,ns=2,ps=instart,data=data.frame(rep(1,3)))
+>
+> mod <- makeDepmix(response=rModels,transition=transition,prior=inMod,ntimes=attr(speed,"ntimes"))
+>
+> ll <- logLik(mod)
+> ll.fb <- logLik(mod,method="fb")
+>
+> logl <- -296.115107102 # see above
+>
+> cat("Test 1: ", all.equal(c(ll),logl,check.att=FALSE), "(loglike of speed data) \n")
+Test 1: TRUE (loglike of speed data)
+>
+>
+> #
+> # model specification made easy
+> #
+>
+> library(depmixS4)
+>
+> resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+> trstart=c(0.899,0.101,0.084,0.916)
+> instart=c(0,1)
+> mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()),respstart=resp,trstart=trstart,instart=instart,prob=T)
+>
+> ll2 <- logLik(mod)
+>
+> cat("Test 1b: ", all.equal(c(ll),c(ll2),check.att=FALSE), "(loglike of speed data) \n")
+Test 1b: TRUE (loglike of speed data)
+>
+> #
+> # TEST 2
+> #
+> # To check the density function for the multinomial responses with a covariate
+> # test a model with a single state, which should be identical to a glm
+> # first fit a model without covariate
+> #
+>
+> invlogit <- function(lp) {exp(lp)/(1+exp(lp))}
+>
+> acc <- glm(corr~1,data=speed,family=binomial)
+>
+> p1 <- invlogit(coef(acc)[1])
+> p0 <- 1-p1
+>
+> mod <- depmix(corr~1,data=speed,nst=1,family=multinomial(),trstart=1,instart=c(1),respstart=c(p0,p1))
+>
+> ll <- logLik(mod)
+>
+> dev <- -2*ll
+>
+> cat("Test 2: ", all.equal(c(dev),acc$deviance),"(loglike of 1-comp glm on acc data) \n")
+Test 2: TRUE (loglike of 1-comp glm on acc data)
+>
+>
+> #
+> # TEST 3
+> #
+> # now add the covariate and compute the loglikelihood
+> #
+>
+> acc <- glm(corr~Pacc,data=speed,family=binomial)
+>
+> p1 <- invlogit(coef(acc)[1])
+> p0 <- 1-p1
+>
+> pstart=c(p0,p1,0,coef(acc)[2])
+>
+> mod <- depmix(corr~Pacc,data=speed,family=multinomial(),trstart=1,instart=1,respst=pstart,nstate=1)
+>
+> ll <- logLik(mod)
+> dev <- -2*ll
+>
+> cat("Test 3: ", all.equal(c(dev),acc$deviance),"(same but now with covariate) \n")
+Test 3: TRUE (same but now with covariate)
+>
+> #
+> # 2-state model with covariate
+> #
+>
+> trstart=c(0.896,0.104,0.084,0.916)
+> trstart=c(trstart[1:2],0,0.01,trstart[3:4],0,0.01)
+> instart=c(0,1)
+> resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+>
+> mod <- depmix(list(rt~1,corr~1),data=speed,family=list(gaussian(),multinomial()),transition=~Pacc,trstart=trstart,instart=instart,respst=resp,nst=2)
+> ll <- logLik(mod)
+>
+> cat("Test 4: ll is now larger than speedll, ie ll is better due to introduction of a covariate \n")
+Test 4: ll is now larger than speedll, ie ll is better due to introduction of a covariate
+> cat("Test 4: ", ll,"\t", logl, "\n")
+Test 4: -296.0333 -296.1151
+> cat("Test 4: ", ll > logl, "\n")
+Test 4: TRUE
+>
Added: trunk/tests/test2getsetpars.R
===================================================================
--- trunk/tests/test2getsetpars.R (rev 0)
+++ trunk/tests/test2getsetpars.R 2008-08-20 14:03:10 UTC (rev 223)
@@ -0,0 +1,20 @@
+#
+# test getpars and setpars
+#
+
+require(depmixS4)
+
+data(speed)
+
+trstart=c(0.896,0.104,0.084,0.916)
+trstart=c(trstart[1:2],0,0.01,trstart[3:4],0,0.01)
+instart=c(0,1)
+resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+
+mod <- depmix(list(rt~1,corr~1),data=speed,family=list(gaussian(),multinomial()),transition=~Pacc,trstart=trstart,instart=instart,respst=resp,nst=2)
+
+mod1 <- setpars(mod,getpars(mod))
+
+cat("Test getpars and setpars: ", all.equal(getpars(mod),getpars(mod1)), "(getpars and setpars) \n")
+
+
Added: trunk/tests/test2getsetpars.Rout.save
===================================================================
--- trunk/tests/test2getsetpars.Rout.save (rev 0)
+++ trunk/tests/test2getsetpars.Rout.save 2008-08-20 14:03:10 UTC (rev 223)
@@ -0,0 +1,41 @@
+
+R version 2.7.1 (2008-06-23)
+Copyright (C) 2008 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+ Natural language support but running in an English locale
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> #
+> # test getpars and setpars
+> #
+>
+> require(depmixS4)
+>
+> data(speed)
+>
+> trstart=c(0.896,0.104,0.084,0.916)
+> trstart=c(trstart[1:2],0,0.01,trstart[3:4],0,0.01)
+> instart=c(0,1)
+> resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+>
+> mod <- depmix(list(rt~1,corr~1),data=speed,family=list(gaussian(),multinomial()),transition=~Pacc,trstart=trstart,instart=instart,respst=resp,nst=2)
+>
+> mod1 <- setpars(mod,getpars(mod))
+>
+> cat("Test getpars and setpars: ", all.equal(getpars(mod),getpars(mod1)), "(getpars and setpars) \n")
+Test getpars and setpars: TRUE (getpars and setpars)
+>
+>
+>
Added: trunk/tests/test3responses.R
===================================================================
--- trunk/tests/test3responses.R (rev 0)
+++ trunk/tests/test3responses.R 2008-08-20 14:03:10 UTC (rev 223)
@@ -0,0 +1,82 @@
+#
+# test response models
+#
+
+require(depmixS4)
+
+
+# gaussian response model
+y <- rnorm(1000)
+mod <- GLMresponse(y~1)
+fm <- fit(mod)
+cat("Test gaussian fit: ", all.equal(getpars(fm),c(mean(y),sd(y)),check=FALSE))
+
+
+# multinomial response model
+y<-sample(1:3,1000,rep=TRUE)
+mod <- GLMresponse(y~1,family=multinomial())
+fm <- fit(mod)
+cat("Test multinomial fit: ", all.equal(getpars(fm),c(table(y)/1000),check=FALSE))
+
+
+x <- sample(0:1,1000,rep=TRUE)
+
+mod <- GLMresponse(sample(1:3,1000,rep=TRUE)~x,family=multinomial(),pstart=c(0.33,0.33,0.33,0,0,1))
+
+mod at y <- simulate(mod)
+
+fit(mod)
+
+colSums(mod at y[which(x==0),])/length(which(x==0))
+
+colSums(mod at y[which(x==1),])/length(which(x==1))
+
+x <- rnorm(1000)
+
+library(boot)
+
+p <- inv.logit(x)
+
+ss <- rbinom(1000,1,p)
+
+mod <- GLMresponse(cbind(ss,1-ss)~x,family=binomial())
+
+fit(mod)
+
+glm(cbind(ss,1-ss)~x, family=binomial)
+
+x <- abs(rnorm(1000,2))
+
+res <- rpois(1000,x)
+
+mod <- GLMresponse(res~x,family=poisson())
+
+fit(mod)
+
+glm(res~x, family=poisson)
+
+x=runif(1000,1,5)
+
+res <- rgamma(1000,x)
+
+## note that gamma needs proper starting values which are not
+## provided by depmixS4 (even with them, this may produce warnings)
+mod <- GLMresponse(res~x,family=Gamma(),pst=c(0.8,1/0.8))
+
+fit(mod)
+
+glm(res~x,family=Gamma)
+
+mn <- c(1,2,3)
+
+sig <- matrix(c(1,.5,0,.5,1,0,0,0,2),3,3)
+
+y <- mvrnorm(1000,mn,sig)
+
+mod <- MVNresponse(y~1)
+
+fit(mod)
+
+colMeans(y)
+
+var(y)
More information about the depmix-commits
mailing list