[Depmix-commits] r644 - in pkg/depmixS4: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 10 12:22:47 CET 2016
Author: ingmarvisser
Date: 2016-02-10 12:22:47 +0100 (Wed, 10 Feb 2016)
New Revision: 644
Modified:
pkg/depmixS4/DESCRIPTION
pkg/depmixS4/R/allGenerics.R
pkg/depmixS4/R/depmix-class.R
pkg/depmixS4/R/depmixfit.R
pkg/depmixS4/R/responseGLM.R
pkg/depmixS4/R/responseGLMMULTINOM.R
pkg/depmixS4/src/fb.h
Log:
=C-header changes
Modified: pkg/depmixS4/DESCRIPTION
===================================================================
--- pkg/depmixS4/DESCRIPTION 2014-12-21 21:21:50 UTC (rev 643)
+++ pkg/depmixS4/DESCRIPTION 2016-02-10 11:22:47 UTC (rev 644)
@@ -7,7 +7,6 @@
Depends: R (>= 3.0.1), nnet, MASS, Rsolnp
Imports: stats, stats4, methods
Suggests: gamlss, gamlss.dist
-Description: Fit latent (hidden) Markov models on mixed categorical and continuous (time series)
- data, otherwise known as dependent mixture models
+Description: Fits latent (hidden) Markov models on mixed categorical and continuous (time series) data, otherwise known as dependent mixture models.
License: GPL (>=2)
URL: http://depmix.r-forge.r-project.org/
Modified: pkg/depmixS4/R/allGenerics.R
===================================================================
--- pkg/depmixS4/R/allGenerics.R 2014-12-21 21:21:50 UTC (rev 643)
+++ pkg/depmixS4/R/allGenerics.R 2016-02-10 11:22:47 UTC (rev 644)
@@ -7,11 +7,11 @@
# version of forward backward routine
.onLoad <- function(lib, pkg) {
- library.dynam("depmixS4", pkg, lib)
+ library.dynam("depmixS4", pkg, lib)
}
.onUnLoad <- function(libpath) {
- library.dynam.unload("depmixS4",libpath)
+ library.dynam.unload("depmixS4",libpath)
}
# Guess what: all generics
Modified: pkg/depmixS4/R/depmix-class.R
===================================================================
--- pkg/depmixS4/R/depmix-class.R 2014-12-21 21:21:50 UTC (rev 643)
+++ pkg/depmixS4/R/depmix-class.R 2016-02-10 11:22:47 UTC (rev 644)
@@ -184,9 +184,9 @@
np <- numeric(object at nresp)
for(j in 1:object at nresp) {
np[j] <- npar(object at response[[1]][[j]])
- pars[[j]] <- matrix(,nr=ns,nc=np[j])
+ pars[[j]] <- matrix(,nrow=ns,ncol=np[j])
}
- allpars <- matrix(,nr=ns,nc=0)
+ allpars <- matrix(,nrow=ns,ncol=0)
nms <- c()
for(j in 1:object at nresp) {
for(i in 1:ns) {
@@ -378,7 +378,7 @@
cat("Transition matrix \n")
pars <- getpars(object)
npprior <- length(getpars(object at prior))
- trm <- matrix(pars[(npprior+1):(ns^2+npprior)],ns,ns,byr=T)
+ trm <- matrix(pars[(npprior+1):(ns^2+npprior)],ns,ns,byrow=T)
if(object at transition[[1]]@family$link == "mlogit") {
trm <- t(apply(trm,1,object at transition[[1]]@family$linkinv,base=object at transition[[1]]@family$base))
}
@@ -436,10 +436,10 @@
for(j in 1:object at nresp) {
#np[j] <- npar(object at response[[1]][[j]]) # this will not always work!
np[j] <- length(allparnames[[j]])
- pars[[j]] <- matrix(,nr=ns,nc=np[j])
+ pars[[j]] <- matrix(,nrow=ns,ncol=np[j])
colnames(pars[[j]]) <- allparnames[[j]]
}
- allpars <- matrix(,nr=ns,nc=0)
+ allpars <- matrix(,nrow=ns,ncol=0)
nms <- c()
for(j in 1:object at nresp) {
for(i in 1:ns) {
Modified: pkg/depmixS4/R/depmixfit.R
===================================================================
--- pkg/depmixS4/R/depmixfit.R 2014-12-21 21:21:50 UTC (rev 643)
+++ pkg/depmixS4/R/depmixfit.R 2016-02-10 11:22:47 UTC (rev 644)
@@ -3,231 +3,227 @@
setMethod("fit",
signature(object="mix"),
function(object,fixed=NULL,equal=NULL,
- conrows=NULL,conrows.upper=NULL,conrows.lower=NULL,
- method=NULL,verbose=TRUE,
- emcontrol=em.control(),
- solnpcntrl=list(rho = 1, outer.iter = 400, inner.iter = 800, delta = 1e-7, tol = 1e-8),
- donlpcntrl=donlp2Control(),
- ...) {
+ conrows=NULL,conrows.upper=NULL,conrows.lower=NULL,
+ method=NULL,verbose=TRUE,
+ emcontrol=em.control(),
+ solnpcntrl=list(rho = 1, outer.iter = 400, inner.iter = 800, delta = 1e-7, tol = 1e-8),
+ donlpcntrl=donlp2Control(),
+ ...) {
- fi <- !is.null(fixed)
- cr <- !is.null(conrows)
- eq <- !is.null(equal)
-
- constr <- any(c(fi,cr,eq))
-
- # when there are constraints donlp/solnp should be used
- # otherwise EM is good
- if(is.null(method)) {
- if(constr) {
- method="rsolnp"
- } else {
- method="EM"
- }
+ fi <- !is.null(fixed)
+ cr <- !is.null(conrows)
+ eq <- !is.null(equal)
+
+ constr <- any(c(fi,cr,eq))
+
+ # when there are constraints donlp/solnp should be used
+ # otherwise EM is good
+ if(is.null(method)) {
+ if(constr) {
+ method="rsolnp"
+ } else {
+ method="EM"
+ }
+ } else {
+ if(method=="EM") {
+ if(constr) {
+ warning("EM not applicable for constrained models; optimization method changed to 'rsolnp'")
+ method="rsolnp"
+ }
+ }
+ }
+
+ if(!(method %in% c("EM","donlp","rsolnp"))) stop("'method' argument invalid; should be one of 'EM', 'rsolnp', 'donlp'.")
+
+ if(method=="EM") {
+ object <- em(object,maxit=emcontrol$maxit,tol=emcontrol$tol,crit=emcontrol$crit,random.start=emcontrol$random.start,classification=emcontrol$classification,verbose=verbose,...)
+ }
+
+ if(method=="donlp"||method=="rsolnp") {
+
+ # check feasibility of starting values
+ if(is.nan(logLik(object))) stop("Initial model infeasible, log likelihood is 'NaN'; please provide better starting values. ")
+
+ # determine which parameters are fixed
+ if(fi) {
+ if(length(fixed)!=npar(object)) stop("'fixed' does not have correct length")
+ } else {
+ if(eq) {
+ if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
+ fixed <- !pa2conr(equal)$free
} else {
- if(method=="EM") {
- if(constr) {
- warning("EM not applicable for constrained models; optimization method changed to 'rsolnp'")
- method="rsolnp"
- }
- }
+ 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]
+
+ constraints <- getConstraints(object)
+
+ lincon=constraints$lincon
+ lin.u=constraints$lin.u
+ lin.l=constraints$lin.l
+ par.u=constraints$par.u
+ par.l=constraints$par.l
+
+ # incorporate equality constraints provided with the fit function, if any
+ if(eq) {
+ 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(cr) {
+ if(ncol(conrows)!=npar(object)) stop("'conrows' does not have the right dimensions")
+ lincon <- rbind(lincon,conrows)
+ if(is.null(conrows.upper)) {
+ 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(is.null(conrows.lower)) {
+ 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,drop=FALSE]
+
+ # remove redundant rows in lincon (all zeroes)
+ allzero <- which(apply(lincon,1,function(y) all(y==0)))
+ if(length(allzero)>0) {
+ lincon <- lincon[-allzero,,drop=FALSE]
+ lin.u <- lin.u[-allzero]
+ lin.l <- lin.l[-allzero]
+ }
+
+ startLogLik <- -logLik(object)*1.01
+
+ # make loglike function that only depends on pars
+ logl <- function(pars) {
+ allpars[!fixed] <- pars
+ object <- setpars(object,allpars)
+ ans = -as.numeric(logLik(object))
+ if(is.na(ans)) ans <- startLogLik
+ if(is.infinite(ans)) ans <- startLogLik
+ ans
+ }
+
+ if(method=="donlp") {
- if(!(method %in% c("EM","donlp","rsolnp"))) stop("'method' argument invalid; should be one of 'EM', 'rsolnp', 'donlp'.")
+ if(!(require(Rdonlp2,quietly=TRUE))) stop("Method 'donlp' requires package 'Rdonlp2'")
+
+ mycontrol <- function(info) {
+ return(TRUE)
+ }
- if(method=="EM") {
- object <- em(object,maxit=emcontrol$maxit,tol=emcontrol$tol,crit=emcontrol$crit,random.start=emcontrol$random.start,classification=emcontrol$classification,verbose=verbose,...)
+ # optimize the parameters
+ result <- donlp2(pars,logl,
+ par.upper=par.u[!fixed],
+ par.lower=par.l[!fixed],
+ A=lincon,
+ lin.upper=lin.u,
+ lin.lower=lin.l,
+ control=donlpcntrl,
+ control.fun=mycontrol,
+ ...
+ )
+
+ if(class(object)=="depmix") object <- as(object,"depmix.fitted") # class(object) <- "depmix.fitted"
+ if(class(object)=="mix") object <- as(object,"mix.fitted") # class(object) <- "mix.fitted"
+
+ # convergence info
+ object at message <- result$message
+
+ # put the result back into the model
+ allpars[!fixed] <- result$par
+ object <- setpars(object,allpars)
+ }
+
+ if(method=="rsolnp") {
+
+ if(!(require(Rsolnp,quietly=TRUE))) stop("Method 'rsolnp' requires package 'Rsolnp'")
+
+ # separate equality and inequality constraints
+ ineq <- which(lin.u!=lin.l)
+ if(length(ineq)>0) lineq <- lincon[-ineq, ,drop=FALSE]
+ else lineq <- lincon
+
+ # returns the evaluated equality constraints
+ if(nrow(lineq)>0) {
+ eqfun <- function(pp) {
+ ans = as.vector(lineq%*%pp)
+ ans
+ }
+ # select the boundary values for the equality constraints
+ if(length(ineq)>0) lineq.bound = lin.l[-ineq]
+ else lineq.bound = lin.l
+ } else {
+ eqfun=NULL
+ lineq.bound=NULL
}
- if(method=="donlp"||method=="rsolnp") {
-
- # check feasibility of starting values
- if(is.nan(logLik(object))) stop("Initial model infeasible, log likelihood is 'NaN'; please provide better starting values. ")
-
- # determine which parameters are fixed
- if(fi) {
- if(length(fixed)!=npar(object)) stop("'fixed' does not have correct length")
- } else {
- if(eq) {
- 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]
-
- constraints <- getConstraints(object)
-
- lincon=constraints$lincon
- lin.u=constraints$lin.u
- lin.l=constraints$lin.l
- par.u=constraints$par.u
- par.l=constraints$par.l
-
- # incorporate equality constraints provided with the fit function, if any
- if(eq) {
- 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(cr) {
- if(ncol(conrows)!=npar(object)) stop("'conrows' does not have the right dimensions")
- lincon <- rbind(lincon,conrows)
- if(is.null(conrows.upper)) {
- 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(is.null(conrows.lower)) {
- 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,drop=FALSE]
-
- # remove redundant rows in lincon (all zeroes)
- allzero <- which(apply(lincon,1,function(y) all(y==0)))
- if(length(allzero)>0) {
- lincon <- lincon[-allzero,,drop=FALSE]
- lin.u <- lin.u[-allzero]
- lin.l <- lin.l[-allzero]
- }
-
- startLogLik <- -logLik(object)*1.01
-
- # make loglike function that only depends on pars
- logl <- function(pars) {
- allpars[!fixed] <- pars
- object <- setpars(object,allpars)
- ans = -as.numeric(logLik(object))
- if(is.na(ans)) ans <- startLogLik
- if(is.infinite(ans)) ans <- startLogLik
- ans
- }
-
- if(method=="donlp") {
-
- reqdon <- require(Rdonlp2,quietly=TRUE)
-
- if(!reqdon) stop("Rdonlp2 not available.")
-
- mycontrol <- function(info) {
- return(TRUE)
- }
-
- # optimize the parameters
- result <- donlp2(pars,logl,
- par.upper=par.u[!fixed],
- par.lower=par.l[!fixed],
- A=lincon,
- lin.upper=lin.u,
- lin.lower=lin.l,
- control=donlpcntrl,
- control.fun=mycontrol,
- ...
- )
-
- if(class(object)=="depmix") object <- as(object,"depmix.fitted") # class(object) <- "depmix.fitted"
- if(class(object)=="mix") object <- as(object,"mix.fitted") # class(object) <- "mix.fitted"
-
- # convergence info
- object at message <- result$message
-
- # put the result back into the model
- allpars[!fixed] <- result$par
- object <- setpars(object,allpars)
- }
-
- if(method=="rsolnp") {
-
- if(!(require(Rsolnp,quietly=TRUE))) stop("Method 'rsolnp' requires package 'Rsolnp'")
-
- # separate equality and inequality constraints
- ineq <- which(lin.u!=lin.l)
- if(length(ineq)>0) lineq <- lincon[-ineq, ,drop=FALSE]
- else lineq <- lincon
-
- # returns the evaluated equality constraints
- if(nrow(lineq)>0) {
- eqfun <- function(pp) {
- ans = as.vector(lineq%*%pp)
- ans
- }
- # select the boundary values for the equality constraints
- if(length(ineq)>0) lineq.bound = lin.l[-ineq]
- else lineq.bound = lin.l
- } else {
- eqfun=NULL
- lineq.bound=NULL
- }
-
- # select the inequality constraints
- if(length(ineq)>0) {
- linineq <- lincon[ineq, ,drop=FALSE]
- ineqLB <- lin.l[ineq]
- ineqUB <- lin.u[ineq]
- } else {
- ineqfun = NULL
- ineqLB=NULL
- ineqUB=NULL
- }
-
- # call to solnp
- res <- solnp(pars,
- logl,
- eqfun = eqfun,
- eqB = lineq.bound,
- ineqfun = ineqfun,
- ineqLB = ineqLB,
- ineqUB = ineqUB,
- LB = par.l[!fixed],
- UB = par.u[!fixed],
- control = solnpcntrl,
- ...
- )
-
- if(class(object)=="depmix") object <- as(object,"depmix.fitted") # class(object) <- "depmix.fitted"
- if(class(object)=="mix") object <- as(object,"mix.fitted") # class(object) <- "mix.fitted"
-
- object at message <- c(res$convergence," (0 is good in Rsolnp, check manual for other values)")
-
- # put the result back into the model
- allpars[!fixed] <- res$pars
- object <- setpars(object,allpars)
-
- }
-
- object at conMat <- linconFull
- object at lin.upper <- lin.u
- object at lin.lower <- lin.l
-
- object at posterior <- viterbi(object)
-
+ # select the inequality constraints
+ if(length(ineq)>0) {
+ linineq <- lincon[ineq, ,drop=FALSE]
+ ineqLB <- lin.l[ineq]
+ ineqUB <- lin.u[ineq]
+ } else {
+ ineqfun = NULL
+ ineqLB=NULL
+ ineqUB=NULL
}
+ # call to solnp
+ res <- solnp(pars,
+ logl,
+ eqfun = eqfun,
+ eqB = lineq.bound,
+ ineqfun = ineqfun,
+ ineqLB = ineqLB,
+ ineqUB = ineqUB,
+ LB = par.l[!fixed],
+ UB = par.u[!fixed],
+ control = solnpcntrl,
+ ...
+ )
+ if(class(object)=="depmix") object <- as(object,"depmix.fitted") # class(object) <- "depmix.fitted"
+ if(class(object)=="mix") object <- as(object,"mix.fitted") # class(object) <- "mix.fitted"
- return(object)
+ object at message <- c(res$convergence," (0 is good in Rsolnp, check manual for other values)")
+
+ # put the result back into the model
+ allpars[!fixed] <- res$pars
+ object <- setpars(object,allpars)
+
+ }
+
+ object at conMat <- linconFull
+ object at lin.upper <- lin.u
+ object at lin.lower <- lin.l
+
+ object at posterior <- viterbi(object)
+
}
+
+ return(object)
+ }
)
@@ -240,7 +236,7 @@
par.l <- rep(-Inf, npar(object))
# make constraint matrix and its upper and lower bounds
- lincon <- matrix(0,nr=0,nc=npar(object))
+ lincon <- matrix(0,nrow=0,ncol=npar(object))
lin.u <- numeric(0)
lin.l <- numeric(0)
Modified: pkg/depmixS4/R/responseGLM.R
===================================================================
--- pkg/depmixS4/R/responseGLM.R 2014-12-21 21:21:50 UTC (rev 643)
+++ pkg/depmixS4/R/responseGLM.R 2016-02-10 11:22:47 UTC (rev 644)
@@ -95,7 +95,7 @@
fixed <- rep(0,ncol(y))
fixed <- c(as.logical(t(fixed)))
constr <- list(
- lin = matrix(1,nr=1,nc=ncol(y)),
+ lin = matrix(1,nrow=1,ncol=ncol(y)),
linup = 1,
linlow = 1,
parup = rep(1,ncol(y)),
Modified: pkg/depmixS4/R/responseGLMMULTINOM.R
===================================================================
--- pkg/depmixS4/R/responseGLMMULTINOM.R 2014-12-21 21:21:50 UTC (rev 643)
+++ pkg/depmixS4/R/responseGLMMULTINOM.R 2016-02-10 11:22:47 UTC (rev 644)
@@ -55,12 +55,12 @@
return(log(rowSums(object at y*predict(object))))
} else {
nr <- nrow(object at y)
- res <- matrix(nr=nr)
+ res <- matrix(nrow=nr)
pr <- predict(object)
# fix this loop!!!! replace with call to apply? or dmultinomial? or other vectorized version?
# possibly use dmultinomial in package mcd2
for(i in 1:nrow(object at y)) {
- res[i,1] <- dmultinom(object at y[i,],pr=pr[i,])
+ res[i,1] <- dmultinom(object at y[i,],prob=pr[i,])
}
return(log(res))
}
@@ -75,12 +75,12 @@
else return(rowSums(object at y*predict(object)))
} else {
nr <- nrow(object at y)
- res <- matrix(nr=nr)
+ res <- matrix(nrow=nr)
pr <- predict(object)
# fix this loop!!!! replace with call to apply? or dmultinomial? or other vectorized version?
# possibly use dmultinomial in package mcd2
for(i in 1:nrow(object at y)) {
- res[i,1] <- dmultinom(object at y[i,],pr=pr[i,])
+ res[i,1] <- dmultinom(object at y[i,],prob=pr[i,])
}
if(log) return(log(res))
else return(res)
Modified: pkg/depmixS4/src/fb.h
===================================================================
--- pkg/depmixS4/src/fb.h 2014-12-21 21:21:50 UTC (rev 643)
+++ pkg/depmixS4/src/fb.h 2016-02-10 11:22:47 UTC (rev 644)
@@ -1,20 +1,14 @@
-
#ifndef FB
#define FB 1
-#include <stdio.h>
-#include <stdlib.h>
-
-extern "C" {
-
-#include <R.h>
-#include <Rmath.h>
+/*
+ * #include <stdio.h>
+ * #include <stdlib.h>
+ */
// compute forward and backward variables, and xi
void forwardbackward(int *hom, int *ns, int *nc, int *nt, int *ntimes, int *bt, int *et,
double *init, double *trdens, double *dens,
double *alpha, double *beta, double *sca, double *xi);
-
-} //end extern "C"
-#endif
\ No newline at end of file
+#endif
More information about the depmix-commits
mailing list