[Depmix-commits] r64 - trunk
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 7 14:12:36 CET 2008
Author: ingmarvisser
Date: 2008-03-07 14:12:36 +0100 (Fri, 07 Mar 2008)
New Revision: 64
Modified:
trunk/depmix-test1balance.R
trunk/depmix-test1balance2.R
trunk/depmix.R
trunk/depmix.fitted.R
trunk/depmixNew-test1.R
trunk/depmixNew-test2.R
trunk/depmixNew-test3.R
trunk/responses.R
Log:
minor changes
Modified: trunk/depmix-test1balance.R
===================================================================
--- trunk/depmix-test1balance.R 2008-03-07 13:11:54 UTC (rev 63)
+++ trunk/depmix-test1balance.R 2008-03-07 13:12:36 UTC (rev 64)
@@ -53,7 +53,7 @@
# 'log Lik.' -1083.036 (df=9)
-source("EM.R")
+source("EM4.R")
# optimize using em
logLik(mod,meth="fb")
Modified: trunk/depmix-test1balance2.R
===================================================================
--- trunk/depmix-test1balance2.R 2008-03-07 13:11:54 UTC (rev 63)
+++ trunk/depmix-test1balance2.R 2008-03-07 13:12:36 UTC (rev 64)
@@ -21,15 +21,20 @@
source("lystig.R")
source("fb.R")
+source("EM.R")
+
# now fit some latent class models
trstart=c(1,0,0,1)
instart=c(0.5,0.5)
# ntimes is added as an argument
+respstart=c(rep(c(0.9,0.1),4),rep(c(0.1,0.9),4))
+
mod <- depmix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
family=list(multinomial(),multinomial(),multinomial(),multinomial()),
- trstart=trstart,instart=instart,ntimes=rep(1,nrow(balance)))
+ trstart=trstart,instart=instart,respstart=respstart,
+ ntimes=rep(1,nrow(balance)))
pars <- getpars(mod)
fixed <- c(1,0,1,1,1,1,rep(c(1,0),8))
Modified: trunk/depmix.R
===================================================================
--- trunk/depmix.R 2008-03-07 13:11:54 UTC (rev 63)
+++ trunk/depmix.R 2008-03-07 13:12:36 UTC (rev 64)
@@ -230,11 +230,11 @@
return(models)
}
-makePriorModel <- function(nstates,ncases,formula=~1,data=NULL,values=NULL) {
+makePriorModel <- function(nstates,ncases,formula=~1,data=NULL,values=NULL, ...) {
# these arguments need to be added at some point FIX ME
base=1
- prob=TRUE
+# prob=TRUE
# initial probabilities model, depending on covariates init(=~1 by default)
if(formula==~1) {
Modified: trunk/depmix.fitted.R
===================================================================
--- trunk/depmix.fitted.R 2008-03-07 13:11:54 UTC (rev 63)
+++ trunk/depmix.fitted.R 2008-03-07 13:12:36 UTC (rev 64)
@@ -30,15 +30,21 @@
setMethod("fit",signature(object="depmix"),
- function(object,w=NULL,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,...) {
+ function(object,w=NULL,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,...) {
# decide on method to be used
# em has to be added here under certain conditions
# in particular when there are linear constraints donlp should be used
- # otherwise donlp is good
- # constraints can not yet be specified however and em is not working properly yet
+ # otherwise EM is good
- method="donlp"
+ # can/does EM deal with fixed constraints??? it should be possible for sure
+ if(is.null(method)) {
+ if(object at stationary&is.null(equal)&is.null(conrows)) {
+ method="EM"
+ } else {
+ method="donlp"
+ }
+ }
# determine which parameters are fixed
if(!is.null(fixed)) {
@@ -53,24 +59,28 @@
}
# set those fixed parameters in the appropriate submodels
object <- setpars(object,fixed,which="fixed")
- # get the full set of parameters
- allpars <- getpars(object)
- # get the reduced set of parameters, ie the ones that will be optimized
- pars <- allpars[!fixed]
-
- # set bounds, if any
- par.u <- rep(+Inf, length(pars))
- par.l <- rep(-Inf, length(pars))
- # make loglike function that only depends on pars
- logl <- function(pars) {
- allpars[!fixed] <- pars
- object <- setpars(object,allpars)
- -logLik(object)
+ if(method=="EM") {
+ object <- em(object,verbose=TRUE)
}
if(method=="donlp") {
+ # get the full set of parameters
+ allpars <- getpars(object)
+ # get the reduced set of parameters, ie the ones that will be optimized
+ pars <- allpars[!fixed]
+ # set bounds, if any
+ par.u <- rep(+Inf, length(pars))
+ par.l <- rep(-Inf, length(pars))
+
+ # make loglike function that only depends on pars
+ logl <- function(pars) {
+ allpars[!fixed] <- pars
+ object <- setpars(object,allpars)
+ -logLik(object)
+ }
+
if(!require(Rdonlp2)) stop("donlp optimization requires the 'Rdonlp2' package")
# make constraint matrix and its upper and lower bounds
@@ -131,7 +141,6 @@
object at conMat <- linconFull
object at message <- result$message
- object at posterior <- matrix(0,1,1) # FIX ME
# put the result back into the model
allpars[!fixed] <- result$par
@@ -139,6 +148,9 @@
}
+ # FIX ME
+ object at posterior <- matrix(0,1,1)
+
return(object)
}
)
@@ -166,6 +178,15 @@
setMethod("show","depmix.fitted",
function(object) {
cat("Convergence info:",object at message,"\n")
+ print(logLik(object))
+ cat("AIC: ", AIC(object),"\n")
+ cat("BIC: ", AIC(object),"\n")
+ }
+)
+
+setMethod("summary","depmix.fitted",
+ function(object) {
+ cat("Convergence info:",object at message,"\n")
cat("Initial state probabilties model \n")
print(object at prior)
cat("\n")
Modified: trunk/depmixNew-test1.R
===================================================================
--- trunk/depmixNew-test1.R 2008-03-07 13:11:54 UTC (rev 63)
+++ trunk/depmixNew-test1.R 2008-03-07 13:12:36 UTC (rev 64)
@@ -14,12 +14,12 @@
# Other tests with optimization of models are moved to depmix-test2.R
#
-setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/code/depmix/trunk/")
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
source("responses.R")
source("depmix.R")
-load("speed.Rda")
+load("data/speed.Rda")
#
# TEST 1: speed data model with optimal parameters, compute the likelihood
@@ -150,9 +150,8 @@
# SPECIFYING MODELS THE EASY WAY
#
-source("depmix.R")
-
mod <- depmix(rt~1,data=speed,nstates=2)
+mod
Modified: trunk/depmixNew-test2.R
===================================================================
--- trunk/depmixNew-test2.R 2008-03-07 13:11:54 UTC (rev 63)
+++ trunk/depmixNew-test2.R 2008-03-07 13:12:36 UTC (rev 64)
@@ -6,19 +6,39 @@
# still works it should return TRUE at every test (or make immediate sense
# otherwise)
-#
-# Test optimization using Rdonlp2
-#
-setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/code/depmix/trunk/")
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
source("responses.R")
source("depmix.R")
source("depmix.fitted.R")
source("llratio.R")
-load("speed.Rda")
+source("EM.R")
+load("data/speed.Rda")
+
+#
+# test model with EM optimization, no covariates
+#
+
+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)
+
+mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()),trstart=trstart)
+# respstart=resp,trstart=trstart,instart=instart)
+
+logLik(mod)
+
+mod1 <- fit(mod)
+ll <- logLik(mod1)
+
+
+#
+# Test optimization using Rdonlp2
+#
+
trstart=c(0.899,0.101,0,0.01,0.084,0.916,0,0)
instart=c(0.5,0.5)
resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
@@ -28,10 +48,6 @@
logLik(mod)
-#
-# fit the model using the general fit function instead
-#
-
mod1 <- fit(mod)
ll <- logLik(mod1)
Modified: trunk/depmixNew-test3.R
===================================================================
--- trunk/depmixNew-test3.R 2008-03-07 13:11:54 UTC (rev 63)
+++ trunk/depmixNew-test3.R 2008-03-07 13:12:36 UTC (rev 64)
@@ -39,37 +39,14 @@
mod1 <- fit(mod,fixed=fixed)
+mod1
+
# 'log Lik.' -1083.036 (df=9)
-logLik(mod1)
-AIC(mod1)
-BIC(mod1)
-
#
# Add age as covariate on class membership
#
-setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
-
-load("data/balance.rda")
-
-source("responses.R")
-source("depmix.R")
-source("depmix.fitted.R")
-
-source("llratio.R")
-source("lystig.R")
-source("fb.R")
-source("EM.R")
-
-# source("responses.R")
-
-# respstart=getpars(mod)[7:22]
-#
-# invlogit <- function(x) {
-# exp(x)/(1+exp(x))
-# }
-# respstart[1:8*2] <- sapply(respstart[1:8*2],invlogit)
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)
Modified: trunk/responses.R
===================================================================
--- trunk/responses.R 2008-03-07 13:11:54 UTC (rev 63)
+++ trunk/responses.R 2008-03-07 13:12:36 UTC (rev 64)
@@ -82,7 +82,7 @@
setMethod("GLMresponse",
signature(formula="formula"),
- function(formula,family,data,pstart=NULL,fixed=NULL) {
+ function(formula,family,data,pstart=NULL,fixed=NULL,prob=TRUE) {
call <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data"), names(mf), 0)
@@ -117,15 +117,19 @@
if(length(pstart)!=npar) stop("length of 'pstart' must be",npar)
if(family$family=="multinomial") {
if(family$link=="identity") parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)])
- else parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)],base=family$base)
+ else {
+# print("ok")
+ if(prob) parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)],base=family$base)
+ else parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
+ }
pstart <- matrix(pstart,ncol(x),byrow=TRUE)
- if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),]
- } else {
+ if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),]
+ } else {
parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)])
}
- if(length(unlist(parameters))>length(parameters$coefficients)) {
+ if(length(unlist(parameters))>length(parameters$coefficients)) {
if(family$family=="gaussian") parameters$sd <- pstart[(length(parameters$coefficients)+1)]
- }
+ }
}
mod <- switch(family$family,
gaussian = new("NORMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar),
@@ -549,6 +553,86 @@
setMethod("fit","transInit",
function(object,w,ntimes) {
+ pars <- object at parameters
+ if(missing(w)) w <- NULL
+ oldfit <- function() {
+ #fit.trMultinom(object,w,ntimes)
+ tol <- 1e-5 # TODO: check global options
+ pars <- object at parameters
+ b <- pars$coefficients
+ base <- object at family$base
+ if(is.matrix(w)) nan <- which(is.na(rowSums(w))) else nan <- which(is.na(w))
+ #vgam(cbind(w[,-base],w[,base]) ~ ) # what is this?
+ y <- as.vector(t(object at family$linkinv(w[-c(nan,ntimes),-base],base=object at family$base)))
+ x <- object at x[-c(nan,ntimes),]
+ if(!is.matrix(x)) x <- matrix(x,ncol=ncol(object at x))
+ nt <- nrow(x)
+ Z <- matrix(ncol=length(b))
+ Z <- vector()
+ for(i in 1:nt) Z <- rbind(Z,t(bdiag(rep(list(x[i,]),ncol(w)-1))))
+ mu <- object at family$linkinv(x%*%b,base=base)
+ mt <- as.numeric(t(mu[,-base]))
+ Dl <- Sigmal <- Wl <- list()
+ converge <- FALSE
+ while(!converge) {
+ b.old <- b
+ for(i in 1:nt) {
+ Dl[[i]] <- object at family$mu.eta(mu[i,-base])
+ Sigmal[[i]] <- object at family$variance(mu[i,-base])
+ Wl[[i]] <- Dl[[i]]%*%solve(Sigmal[[i]])%*%t(Dl[[i]]) # TODO:
+ }
+ Sigma <- bdiag(Sigmal)
+ D <- bdiag(Dl)
+ W <- bdiag(Wl)
+
+ b[,-base] <- as.numeric(b[,-base]) + solve(t(Z)%*%W%*%Z)%*%(t(Z)%*%D%*%solve(Sigma)%*%(y-mt))
+ if(abs(sum(b-b.old)) < tol) converge <- TRUE
+ mu <- object at family$linkinv(x%*%b,base=base)
+ mt <- as.numeric(t(mu[,-base]))
+ }
+ pars$coefficients <- t(b) # TODO: setpars gets matrix in wrong order!!! Fix this in setpars.
+ pars
+ }
+
+ vglmfit <- function() {
+ base <- object at family$base
+ w <- cbind(w[,-base],w[,base])
+ x <- slot(object,"x")
+ fam <- slot(object,"family")
+ fit <- vglm(w~x,fam)
+ pars$coefficients[,-base] <- t(slot(fit,coefficients)) # TODO: setpars gets matrix in wrong order!!! Fix this in setpars.
+ pars
+ }
+
+ nnetfit <- function() {
+ require(nnet)
+ pars <- object at parameters
+ base <- object at family$base # delete me
+ #y <- object at y[,-base]
+ y <- object at y
+ x <- object at x
+ if(is.matrix(y)) na <- unlist(apply(y,2,function(x) which(is.na(x)))) else na <- which(is.na(y))
+ if(is.matrix(x)) na <- c(na,unlist(apply(x,2,function(x) which(is.na(x))))) else na <- c(na,which(is.na(x)))
+ if(!is.null(w)) na <- c(na,which(is.na(w)))
+ y <- as.matrix(y)
+ x <- as.matrix(x)
+ na <- unique(na)
+ x <- x[-na,]
+ y <- y[-na,]
+ y <- round(y) # delete me
+ if(!is.null(w)) w <- w[-na]
+ #mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
+ #mask[,base] <- 0
+ if(!is.null(w)) fit <- multinom(y~x-1,weights=w,trace=FALSE) else fit <- multinom(y~x-1,weights=w,trace=FALSE)
+ ids <- vector(,length=ncol(y))
+ ids[base] <- 1
+ ids[-base] <- 2:ncol(y)
+ pars$coefficients <- t(matrix(fit$wts,ncol=ncol(y))[-1,ids])
+ object <- setpars(object,unlist(pars))
+ #object
+ pars
+ }
+
require(nnet)
pars <- object at parameters
base <- object at family$base # delete me
@@ -565,8 +649,8 @@
x <- x[-na,]
y <- y[-na,]
#y <- round(y) # delete me
- if(!is.null(w)) w <- w[-na]
- }
+ if(!is.null(w)) w <- w[-na]
+ }
#mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
#mask[,base] <- 0
if(!is.null(w)) fit <- multinom(y~x-1,weights=w,trace=FALSE) else fit <- multinom(y~x-1,trace=FALSE)
@@ -576,6 +660,6 @@
pars$coefficients <- t(matrix(fit$wts,ncol=ncol(y))[-1,ids]) # why do we need to transpose?
object <- setpars(object,unlist(pars))
object
+ }
+)
- }
-)
\ No newline at end of file
More information about the depmix-commits
mailing list