SCM

SCM Repository

[depmix] Diff of /trunk/R/depmix.fitted.R
ViewVC logotype

Diff of /trunk/R/depmix.fitted.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 20, Sun Mar 2 09:03:06 2008 UTC revision 64, Fri Mar 7 13:12:36 2008 UTC
# Line 1  Line 1 
  # # Ingmar Visser, 29-2-2008 # # # DEPMIX.FITTED CLASS # # this has information from fitting a depmix model # 1) whether fitting was succesful: convergence etc # 2) some gof information: logLik, AIC, BIC etc: these are provided through methods on the depmix object # 3) constraints that were fitted on the parameters # 4) posterior probabilities or viterbi sequence 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        ),        contains="depmix" ) # Which methods are needed for this class? # 1) construction method: the fit function which is called on a depmix # class object, this is where the real work is done: done!!! # 2) print method: done!!!! # 3) summary method setMethod("fit",signature(object="depmix"),        function(object,w=NULL,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,...) {                                # 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                                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")                # 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=="donlp") {                                                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@conMat <- linconFull                        object@message <- result$message                        object@posterior <- matrix(0,1,1) # FIX ME                                                # put the result back into the model                        allpars[!fixed] <- result$par                        object <- setpars(object,allpars)                                        }                                return(object)        } ) # convert conpat vector to rows of constraint matrix pa2conr <- function(x) {        fix=as.logical(x)        x=replace(x,list=which(x==1),0)        un=setdiff(unique(x),0)        y=matrix(0,0,length(x))        for(i in un) {                z=which(x==i)                for(j in 2:length(z)) {                        k=rep(0,length(x))                        k[z[1]]=1                        k[z[j]]=-1                        y=rbind(y,k)                }        }        pa = list(free=fix,conr=y)        return(pa) } setMethod("show","depmix.fitted",        function(object) {                cat("Convergence info:",object@message,"\n")                cat("Initial state probabilties model \n")                print(object@prior)                cat("\n")                for(i in 1:object@nstates) {                        cat("Transition model for state (component)", i,"\n")                        print(object@transition[[i]])                        cat("\n")                }                cat("\n")                for(i in 1:object@nstates) {                        cat("Response model(s) for state", i,"\n\n")                        for(j in 1:object@nresp) {                                cat("Response model for response",j,"\n")                                print(object@response[[i]][[j]])                                cat("\n")                        }                        cat("\n")                }        } )  
1     # # Ingmar Visser, 29-2-2008 # # # DEPMIX.FITTED CLASS # # this has information from fitting a depmix model # 1) whether fitting was succesful: convergence etc # 2) some gof information: logLik, AIC, BIC etc: these are provided through methods on the depmix object # 3) constraints that were fitted on the parameters # 4) posterior probabilities or viterbi sequence 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        ),        contains="depmix" ) # Which methods are needed for this class? # 1) construction method: the fit function which is called on a depmix # class object, this is where the real work is done: done!!! # 2) print method: done!!!! # 3) summary method setMethod("fit",signature(object="depmix"),        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 EM is good                                # can/does EM deal with fixed constraints??? it should be possible for sure                if(is.null(method)) {                        if(object@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@conMat <- linconFull                        object@message <- result$message                                                # put the result back into the model                        allpars[!fixed] <- result$par                        object <- setpars(object,allpars)                                        }                                # FIX ME                object@posterior <- matrix(0,1,1)                                return(object)        } ) # convert conpat vector to rows of constraint matrix pa2conr <- function(x) {        fix=as.logical(x)        x=replace(x,list=which(x==1),0)        un=setdiff(unique(x),0)        y=matrix(0,0,length(x))        for(i in un) {                z=which(x==i)                for(j in 2:length(z)) {                        k=rep(0,length(x))                        k[z[1]]=1                        k[z[j]]=-1                        y=rbind(y,k)                }        }        pa = list(free=fix,conr=y)        return(pa) } setMethod("show","depmix.fitted",        function(object) {                cat("Convergence info:",object@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@message,"\n")                cat("Initial state probabilties model \n")                print(object@prior)                cat("\n")                for(i in 1:object@nstates) {                        cat("Transition model for state (component)", i,"\n")                        print(object@transition[[i]])                        cat("\n")                }                cat("\n")                for(i in 1:object@nstates) {                        cat("Response model(s) for state", i,"\n\n")                        for(j in 1:object@nresp) {                                cat("Response model for response",j,"\n")                                print(object@response[[i]][[j]])                                cat("\n")                        }                        cat("\n")                }        } )

Legend:
Removed from v.20  
changed lines
  Added in v.64

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge