[Depmix-commits] r126 - pkg trunk trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 20 09:31:27 CEST 2008
Author: ingmarvisser
Date: 2008-05-20 09:31:27 +0200 (Tue, 20 May 2008)
New Revision: 126
Added:
trunk/R/depmixAIC.R
trunk/R/depmixBIC.R
trunk/R/depmixfit-class.R
trunk/R/depmixfit.R
trunk/R/transInit.R
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
trunk/NAMESPACE
Log:
Put most functions into separate R files (part 2)
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2008-05-19 20:29:57 UTC (rev 125)
+++ pkg/DESCRIPTION 2008-05-20 07:31:27 UTC (rev 126)
@@ -1,6 +1,6 @@
Package: depmixS4
-Version: 0.1-0
-Date: 2008-03-23
+Version: 0.1-1
+Date: 2008-05-19
Title: Dependent Mixture Models
Author: Ingmar Visser <i.visser at uva.nl>, Maarten Speekenbrink <m.speekenbrink at ucl.ac.uk>
Maintainer: Ingmar Visser <i.visser at uva.nl>
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2008-05-19 20:29:57 UTC (rev 125)
+++ pkg/NAMESPACE 2008-05-20 07:31:27 UTC (rev 126)
@@ -1,7 +1,6 @@
import(methods)
importFrom(stats, predict)
-importFrom(stats4, AIC, BIC, logLik, plot, summary)
export(
makeDepmix,
Modified: trunk/NAMESPACE
===================================================================
--- trunk/NAMESPACE 2008-05-19 20:29:57 UTC (rev 125)
+++ trunk/NAMESPACE 2008-05-20 07:31:27 UTC (rev 126)
@@ -1,7 +1,6 @@
import(methods)
importFrom(stats, predict)
-importFrom(stats4, AIC, BIC, logLik, plot, summary)
export(
makeDepmix,
Added: trunk/R/depmixAIC.R
===================================================================
--- trunk/R/depmixAIC.R (rev 0)
+++ trunk/R/depmixAIC.R 2008-05-20 07:31:27 UTC (rev 126)
@@ -0,0 +1,6 @@
+# depends on logLik and freepars
+setMethod("AIC", signature(object="depmix"),
+ function(object, ..., k=2){
+ c(-2 * logLik(object) + freepars(object) * k)
+ }
+)
Added: trunk/R/depmixBIC.R
===================================================================
--- trunk/R/depmixBIC.R (rev 0)
+++ trunk/R/depmixBIC.R 2008-05-20 07:31:27 UTC (rev 126)
@@ -0,0 +1,6 @@
+# depends on logLik, freepars and nobs
+setMethod("BIC", signature(object="depmix"),
+ function(object, ...){
+ c(-2 * logLik(object) + freepars(object) * log(nobs(object)))
+ }
+)
Added: trunk/R/depmixfit-class.R
===================================================================
--- trunk/R/depmixfit-class.R (rev 0)
+++ trunk/R/depmixfit-class.R 2008-05-20 07:31:27 UTC (rev 126)
@@ -0,0 +1,58 @@
+
+#
+# Ingmar Visser, 23-3-2008
+#
+
+#
+# DEPMIX.FITTED CLASS
+#
+
+setClass("depmix.fitted",
+ representation(message="character", # convergence information
+ conMat="matrix", # constraint matrix on the parameters for general linear constraints
+ posterior="data.frame" # posterior probabilities for the states
+ ),
+ contains="depmix"
+)
+
+# accessor function
+
+setMethod("posterior","depmix.fitted",
+ function(object) {
+ return(object at posterior)
+ }
+)
+
+setMethod("show","depmix.fitted",
+ function(object) {
+ cat("Convergence info:",object at message,"\n")
+ print(logLik(object))
+ cat("AIC: ", AIC(object),"\n")
+ cat("BIC: ", BIC(object),"\n")
+ }
+)
+
+setMethod("summary","depmix.fitted",
+ function(object) {
+ cat("Initial state probabilties model \n")
+ print(object at prior)
+ cat("\n")
+ for(i in 1:object at nstates) {
+ cat("Transition model for state (component)", i,"\n")
+ print(object at transition[[i]])
+ cat("\n")
+ }
+ cat("\n")
+ for(i in 1:object at nstates) {
+ cat("Response model(s) for state", i,"\n\n")
+ for(j in 1:object at nresp) {
+ cat("Response model for response",j,"\n")
+ print(object at response[[i]][[j]])
+ cat("\n")
+ }
+ cat("\n")
+ }
+ }
+)
+
+
Added: trunk/R/depmixfit.R
===================================================================
--- trunk/R/depmixfit.R (rev 0)
+++ trunk/R/depmixfit.R 2008-05-20 07:31:27 UTC (rev 126)
@@ -0,0 +1,125 @@
+
+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
+
+ # 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)) {
+ if(length(fixed)!=npar(object)) stop("'fixed' does not have correct length")
+ } else {
+ if(!is.null(equal)) {
+ if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
+ fixed <- !pa2conr(equal)$free
+ } else {
+ fixed <- getpars(object,"fixed")
+ }
+ }
+ # set those fixed parameters in the appropriate submodels
+ object <- setpars(object,fixed,which="fixed")
+
+ 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
+ lincon <- matrix(0,nr=0,nc=npar(object))
+ lin.u <- numeric(0)
+ lin.l <- numeric(0)
+
+ # incorporate equality constraints, if any
+ if(!is.null(equal)) {
+ if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
+ equal <- pa2conr(equal)$conr
+ lincon <- rbind(lincon,equal)
+ lin.u <- c(lin.u,rep(0,nrow(equal)))
+ lin.l <- c(lin.l,rep(0,nrow(equal)))
+ }
+
+ # incorporate general linear constraints, if any
+ if(!is.null(conrows)) {
+ if(ncol(conrows)!=npar(object)) stop("'conrows' does not have the right dimensions")
+ lincon <- rbind(lincon,conrows)
+ if(conrows.upper==0) {
+ lin.u <- c(lin.u,rep(0,nrow(conrows)))
+ } else {
+ if(length(conrows.upper)!=nrow(conrows)) stop("'conrows.upper does not have correct length")
+ lin.u <- c(lin.u,conrows.upper)
+ }
+ if(conrows.lower==0) {
+ lin.l <- c(lin.l,rep(0,nrow(conrows)))
+ } else {
+ if(length(conrows.lower)!=nrow(conrows)) stop("'conrows.lower does not have correct length")
+ lin.l <- c(lin.l,conrows.lower)
+ }
+ }
+
+ # select only those columns of the constraint matrix that correspond to non-fixed parameters
+ linconFull <- lincon
+ lincon <- lincon[,!fixed]
+
+ # set donlp2 control parameters
+ cntrl <- donlp2.control(hessian=FALSE,difftype=2,report=TRUE)
+
+ mycontrol <- function(info) {
+ return(TRUE)
+ }
+
+ # optimize the parameters
+ result <- donlp2(pars,logl,
+ par.upper=par.u,
+ par.lower=par.l,
+ A=lincon,
+ lin.upper=lin.u,
+ lin.lower=lin.l,
+ control=cntrl,
+ control.fun=mycontrol,
+ ...
+ )
+
+ class(object) <- "depmix.fitted"
+
+ object at conMat <- linconFull
+ object at message <- result$message
+
+ # put the result back into the model
+ allpars[!fixed] <- result$par
+ object <- setpars(object,allpars)
+
+ }
+
+ object at posterior <- viterbi2(object)
+
+ return(object)
+ }
+)
\ No newline at end of file
Added: trunk/R/transInit.R
===================================================================
--- trunk/R/transInit.R (rev 0)
+++ trunk/R/transInit.R 2008-05-20 07:31:27 UTC (rev 126)
@@ -0,0 +1,189 @@
+#
+# for the transition models and the prior (y is missing, ie there is no
+# response, and nstates must be provided as the number of categories
+# neccessary in the mulinomial model)
+#
+
+setClass("transInit",contains="GLMresponse")
+
+# FIX ME: data is a necessary argument to determine the dimension of x, even when there
+# are no covariates (and there are by definition no responses ...)
+setMethod("transInit",
+ signature(formula="formula"),
+ function(formula,nstates,data=NULL,family=multinomial(),pstart=NULL,fixed=NULL,prob=TRUE, ...) {
+ call <- match.call()
+ mf <- match.call(expand.dots = FALSE)
+ m <- match(c("formula", "data"), names(mf), 0)
+ mf <- mf[c(1, m)]
+ mf$drop.unused.levels <- TRUE
+ mf[[1]] <- as.name("model.frame")
+ mf <- eval(mf, parent.frame())
+ x <- model.matrix(attr(mf, "terms"),mf)
+ y <- matrix(1,ncol=1) # y is not needed in the transition and init models
+ parameters <- list()
+ if(is.null(nstates)) stop("'nstates' must be provided in call to trinModel")
+ if(family$family=="multinomial") {
+ parameters$coefficients <- matrix(0,ncol=nstates,nrow=ncol(x))
+ if(is.null(fixed)) {
+ fixed <- parameters$coefficients
+ fixed[,family$base] <- 1
+ fixed <- c(as.logical(t(fixed)))
+ }
+ }
+ npar <- length(unlist(parameters))
+ if(is.null(fixed)) fixed <- rep(0,npar)
+ if(!is.null(pstart)) {
+ if(length(pstart)!=npar) stop("length of 'pstart' must be ",npar)
+ if(family$family=="multinomial") {
+ if(prob) {
+ 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 {
+ 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(family$link=="identity") parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)])
+ else parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)],base=family$base)
+ }
+ }
+ mod <- switch(family$family,
+ multinomial = new("transInit",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar),
+ new("transInit",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar)
+ )
+ mod
+ }
+)
+
+setMethod("logDens","transInit",
+ function(object) {
+ log(predict(object))
+ }
+)
+
+setMethod("dens","transInit",
+ function(object,log=FALSE) {
+ if(log) log(predict(object))
+ else predict(object)
+ }
+)
+
+setMethod("predict","transInit",
+ function(object) {
+ object at family$linkinv(object at x%*%object at parameters$coefficients,base=object at family$base)
+ }
+)
+
+setMethod("fit","transInit",
+ function(object,w,ntimes) {
+ pars <- object at parameters
+ if(missing(w)) w <- NULL
+ oldfit <- function() {
+ 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() {
+ 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
+ }
+
+ 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)
+ if(length(na)>0) {
+ 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,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]) # why do we need to transpose?
+ object <- setpars(object,unlist(pars))
+ object
+ }
+)
+
More information about the depmix-commits
mailing list