[Depmix-commits] r352 - in trunk: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 22 14:28:34 CET 2010
Author: ingmarvisser
Date: 2010-02-22 14:28:34 +0100 (Mon, 22 Feb 2010)
New Revision: 352
Modified:
trunk/R/depmixfit.R
trunk/depmixNew-test4.R
trunk/depmixNew-test5.R
trunk/depmixNew-test6.R
Log:
Added Rsolnp support including an example in ...test6.R
Modified: trunk/R/depmixfit.R
===================================================================
--- trunk/R/depmixfit.R 2010-02-22 13:25:47 UTC (rev 351)
+++ trunk/R/depmixfit.R 2010-02-22 13:28:34 UTC (rev 352)
@@ -1,131 +1,208 @@
-
-setMethod("fit",
- signature(object="mix"),
- function(object,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,...) {
-
- # when there are 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(is.null(equal)&is.null(conrows)&is.null(fixed)) {
- 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(is.nan(logLik(object))) stop("Initial model infeasible, log likelihood is NaN; please provide better starting values. ")
-
- 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(any(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(any(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,drop=FALSE]
-
- # 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,
- ...
- )
-
- if(class(object)=="depmix") class(object) <- "depmix.fitted"
- if(class(object)=="mix") class(object) <- "mix.fitted"
-
- object at conMat <- linconFull
- object at message <- result$message
- object at lin.upper <- lin.u
- object at lin.lower <- lin.l
-
- # put the result back into the model
- allpars[!fixed] <- result$par
- object <- setpars(object,allpars)
-
- }
-
- object at posterior <- viterbi(object)
-
- return(object)
- }
+
+setMethod("fit",
+ signature(object="mix"),
+ function(object,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,...) {
+
+ 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="donlp"
+ } else {
+ method="EM"
+ }
+ } else {
+ if(method=="EM") {
+ if(constr) {
+ warning("EM not applicable for constrained models; optimization method changed to 'donlp'")
+ method="donlp"
+ }
+ }
+ }
+
+
+ # check feasibility of starting values
+ if(is.nan(logLik(object))) stop("Initial model infeasible, log likelihood is NaN; please provide better starting values. ")
+
+ if(method=="EM") {
+ object <- em(object,verbose=TRUE,...)
+ }
+
+ if(method=="donlp"||method=="rsolnp") {
+
+ if(method=="donlp"&!require(Rdonlp2)) method="rsolnp"
+
+ if(method=="rsolnp"&!(require(Rsolnp))) stop("Optimization requires either 'Rdonlp2' or 'Rsolnp'")
+
+ # 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]
+
+ # set bounds, if any (should add bounds for eg sd parameters at some point ...)
+ 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)
+ ans = -as.numeric(logLik(object))
+ if(is.na(ans)) ans = 100000
+ ans
+ }
+
+ # 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(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(any(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(any(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,drop=FALSE]
+
+ if(method=="donlp") {
+ # 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,
+ ...
+ )
+
+ # 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") {
+
+ # 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,
+ eqgrad =NULL,
+ ineqfun = ineqfun,
+ ineqLB = ineqLB,
+ ineqUB = ineqUB,
+ ineqgrad = NULL,
+ LB = NULL,
+ UB = NULL,
+ control = list(trace = 1)
+ )
+
+ 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)
+
+ }
+
+ if(class(object)=="depmix") class(object) <- "depmix.fitted"
+ if(class(object)=="mix") class(object) <- "mix.fitted"
+
+ object at conMat <- linconFull
+ object at lin.upper <- lin.u
+ object at lin.lower <- lin.l
+
+ }
+
+ object at posterior <- viterbi(object)
+
+ return(object)
+ }
)
\ No newline at end of file
Modified: trunk/depmixNew-test4.R
===================================================================
--- trunk/depmixNew-test4.R 2010-02-22 13:25:47 UTC (rev 351)
+++ trunk/depmixNew-test4.R 2010-02-22 13:28:34 UTC (rev 352)
@@ -11,7 +11,7 @@
# library(depmixS4)
-# setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
+# setwd("~/Documents/projects/depmixProject/codesvn/depmix/trunk/")
#
# optimization speed profile: case 1: latent class data
Modified: trunk/depmixNew-test5.R
===================================================================
--- trunk/depmixNew-test5.R 2010-02-22 13:25:47 UTC (rev 351)
+++ trunk/depmixNew-test5.R 2010-02-22 13:28:34 UTC (rev 352)
@@ -85,96 +85,3 @@
-
-#
-# multivariate normal mixture models
-#
-
-
-
-library(depmixS4)
-# use function xpnd and vech from MCMCpack to convert from lower.tri to square matrix and back
-
-# multivariate normal response model
-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~rnorm(1000))
-
-y <- simulate(mod)
-
-head(dens(mod,log=T))
-
-head(predict(mod))
-
-mod <- fit(mod)
-colMeans(y)
-var(y)
-
-mod
-
-npar(mod)
-
-require(MASS)
-
-m1 <- c(0,1)
-sd1 <- matrix(c(1,0.7,.7,1),2,2)
-
-m2 <- c(1,0)
-sd2 <- matrix(c(2,.1,.1,1),2,2)
-
-y1 <- mvrnorm(50,m1,sd1)
-y2 <- mvrnorm(50,m2,sd2)
-
-y <- rbind(y1,y2)
-
-m1 <- MVNresponse(y~1,pst=c(0,.1,1,0.1,1))
-
-m2 <- MVNresponse(y~1)
-
-m1
-
-m1 at parameters
-
-m2
-
-m2 at parameters
-
-rModels <- list(
- list(
- MVNresponse(y~1)
- ),
- list(
- MVNresponse(y~1)
- )
-)
-
-trstart=c(0.9,0.1,0.1,0.9)
-
-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=runif(2)
-inMod <- transInit(~1,ns=2,ps=instart,data=data.frame(1))
-
-mod <- makeDepmix(response=rModels,transition=transition,prior=inMod)
-
-logLik(mod)
-
-
-fm <- fit(mod)
-fmd <- fit(mod,meth="donlp")
-
-pem <- getpars(fm)[7:16]
-pdon <- getpars(fmd)[7:16]
-
-all.equal(pem,pdon)
-
-fm <- simulate(fm)
-
-fm <- fit(fm)
-
-fm
-
-summary(fm)
Modified: trunk/depmixNew-test6.R
===================================================================
--- trunk/depmixNew-test6.R 2010-02-22 13:25:47 UTC (rev 351)
+++ trunk/depmixNew-test6.R 2010-02-22 13:28:34 UTC (rev 352)
@@ -7,30 +7,6 @@
setwd("~/Documents/projects/depmixProject/codesvn/depmix/trunk/")
-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)
-}
-
-
-#
-# Rsolnp optimization
-#
-
require(depmixS4)
data(speed)
@@ -45,10 +21,6 @@
logLik(mod1)
-# mod1 <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,nstates=2,family=list(gaussian(),multinomial()),
-# trstart=trst,respst=c(5,.5,.5,.5,6,.5,.1,.9))
-# logLik(mod1)
-
# fit the model
fmod1 <- fit(mod1)
fmod1 # to see the logLik and optimization information
@@ -56,92 +28,38 @@
# to see the parameters
summary(fmod1)
-
-require(Rsolnp)
-
allpars <- getpars(fmod1)
allpars[2]=Inf # this means the process will always start in state 2
allpars[14]=0 # the corr parameters in state 1 are now both 0, corresponding the 0.5 prob
-
-lowb <- rep(-Inf,length(allpars))
-lowb[c(12,16)] <- 0
-
-
-fixed <- getpars(fmod1,"fixed")
-fixed[2] <- TRUE
-fixed[14] <- TRUE
-
-pars <- allpars[!fixed]
-lowb <- lowb[!fixed]
-
-set.seed(2)
-
-pars <- pars+runif(9,0,.1)
-allpars[!fixed] <- pars
-stmod <- setpars(fmod1,allpars)
-
-fm1Rdon <- fit(stmod,fixed=fixed)
-
-
-
-# define loglike function
-ll <- function(pars) {
- allpars[!fixed] <- pars
- mod1 <- setpars(mod1,allpars)
- ans = -as.numeric(logLik(mod1))
- if(is.na(ans)) ans = 100000
- ans
-}
-
-ll(pars)
-
-# use Rsolnp for unconstrained model
-res <- solnp(pars, ll, control = list(trace = 1),LB=lowb)
-
-res$pars
-
-fpars <- allpars
-fpars[!fixed] <- res$pars
-
-fmRsol <- setpars(mod1,fpars)
-
-
-
-
allpars[c(4,8)] <- -4
allpars[c(6,10)] <- 10
+stmod <- setpars(fmod1,allpars)
+
conpat <- c(0,0,rep(c(0,1),4),1,1,0,0,1,1,0,1)
# constrain the beta's on the transition parameters to be equal
conpat[4] <- conpat[8] <- 2
conpat[6] <- conpat[10] <- 3
-eqA=pa2conr(conpat)$conr
-eqA.bound=c(0,0)
+fm1sol <- fit(stmod,equal=conpat,method="rsolnp")
+fm1don <- fit(stmod,equal=conpat,method="donlp")
-eqA <- eqA[,!fixed]
-lowb <- rep(-Inf,length(allpars))
-lowb[c(12,16)] <- 0
+fm1sol
+fm1don
-pars <- allpars[!fixed]
-lowb <- lowb[!fixed]
+summary(fm1sol)
+summary(fm1don)
-eqfun <- function(x) {
- ans = as.vector(eqA %*% x)
- ans
-}
+getpars(fm1sol)
+getpars(fm1don)
-eqfun(pars)
-ll(pars)
-res <- solnp(pars, ll, eqfun = eqfun, eqB = eqA.bound,LB=lowb)
-fm1 <- fit(mod1,conpat=conpat)
@@ -154,8 +72,3 @@
-
-
-
-
-
More information about the depmix-commits
mailing list