[Depmix-commits] r444 - in pkg/depmixS4: . R inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 20 13:01:22 CEST 2010
Author: ingmarvisser
Date: 2010-10-20 13:01:22 +0200 (Wed, 20 Oct 2010)
New Revision: 444
Added:
pkg/depmixS4/R/em.control.R
pkg/depmixS4/man/em.control.Rd
Modified:
pkg/depmixS4/NAMESPACE
pkg/depmixS4/NEWS
pkg/depmixS4/R/EM.R
pkg/depmixS4/R/depmixfit.R
pkg/depmixS4/inst/doc/depmixS4.Rnw
pkg/depmixS4/man/depmix.fit.Rd
pkg/depmixS4/man/makeDepmix.Rd
Log:
Added function em.control that returns the list of em control parameters. This makes future additions to this list easier and gives less clutter in the function definition of fit. The default for EM is now to use random start values.
Modified: pkg/depmixS4/NAMESPACE
===================================================================
--- pkg/depmixS4/NAMESPACE 2010-10-19 15:08:36 UTC (rev 443)
+++ pkg/depmixS4/NAMESPACE 2010-10-20 11:01:22 UTC (rev 444)
@@ -11,7 +11,8 @@
MVNresponse,
llratio,
multinomial,
- em,
+ em,
+ em.control,
viterbi,
mlogit
)
Modified: pkg/depmixS4/NEWS
===================================================================
--- pkg/depmixS4/NEWS 2010-10-19 15:08:36 UTC (rev 443)
+++ pkg/depmixS4/NEWS 2010-10-20 11:01:22 UTC (rev 444)
@@ -1,4 +1,20 @@
+Changes in depmixS4 version 1.0-1
+
+ o minor changes in documentation to conform to R 2.12.0 standards.
+
+ o fixed a bug concerning random start values (the argument to specify
+ this was not passed to the EM algorithm ...).
+
+ o changed the emcontrol argument to the fit function; it now calls
+ a function em.control which returns the list of control parameters, which
+ now also includes maxit, the max number of iterations of the EM algorithm.
+ This makes future additions to EM control parameters easier.
+
+ o The use of random parameter initialization is now the default when using EM
+ to fit models.
+
+
Changes in depmixS4 version 1.0-0
o added a vignette to the package and upped the version number 1.0-0 to
Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R 2010-10-19 15:08:36 UTC (rev 443)
+++ pkg/depmixS4/R/EM.R 2010-10-20 11:01:22 UTC (rev 444)
@@ -15,12 +15,10 @@
}
# em for lca and mixture models
-em.mix <- function(object,maxit=100,tol=1e-8,crit=c("relative","absolute"),random.start=FALSE,verbose=FALSE,...) {
+em.mix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,...) {
if(!is(object,"mix")) stop("object is not of class 'mix'")
-
- crit <- match.arg(crit)
-
+
ns <- nstates(object)
ntimes <- ntimes(object)
lt <- length(ntimes)
@@ -30,12 +28,12 @@
converge <- FALSE
j <- 0
- # compute responsibilities
+ # compute response probabilities
B <- apply(object at dens,c(1,3),prod)
gamma <- object at init*B
LL <- sum(log(rowSums(gamma)))
# normalize
- gamma <- gamma/rowSums(gamma)
+ gamma <- gamma/rowSums(gamma)
if(random.start) {
nr <- sum(ntimes(object))
@@ -110,10 +108,9 @@
}
# em for hidden markov models
-em.depmix <- function(object,maxit=100,tol=1e-8,crit=c("relative","absolute"),random.start=FALSE,verbose=FALSE,...) {
+em.depmix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,...) {
- if(!is(object,"depmix")) stop("object is not of class '(dep)mix'")
- crit <- match.arg(crit)
+ if(!is(object,"depmix")) stop("object is not of class 'depmix'")
ns <- nstates(object)
@@ -125,10 +122,6 @@
converge <- FALSE
j <- 0
- # A <- object at trDens
- # B <- object at dens
- # init <- object at init
-
# initial expectation
fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
LL <- fbo$logLike
@@ -144,9 +137,9 @@
# maximization
- # should become object at prior <- fit(object at prior)
+ # should become object at prior <- fit(object at prior, gamma)
object at prior@y <- fbo$gamma[bt,,drop=FALSE]
- object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
+ object at prior <- fit(object at prior, w=NULL, ntimes=NULL)
object at init <- dens(object at prior)
trm <- matrix(0,ns,ns)
@@ -158,7 +151,8 @@
for(k in 1:ns) {
trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
}
- # FIX THIS; it will only work with specific trinModels??
+ # FIX THIS; it will only work with specific trinModels
+ # should become object at transition = fit(object at transition, xi, gamma)
object at transition[[i]]@parameters$coefficients <- switch(object at transition[[i]]@family$link,
identity = object at transition[[i]]@family$linkfun(trm[i,]),
mlogit = object at transition[[i]]@family$linkfun(trm[i,],base=object at transition[[i]]@family$base),
Modified: pkg/depmixS4/R/depmixfit.R
===================================================================
--- pkg/depmixS4/R/depmixfit.R 2010-10-19 15:08:36 UTC (rev 443)
+++ pkg/depmixS4/R/depmixfit.R 2010-10-20 11:01:22 UTC (rev 444)
@@ -1,7 +1,8 @@
+
setMethod("fit",
signature(object="mix"),
- function(object,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,em.control=list(tol=1e-8,crit=c("relative","absolute"),random.start=FALSE),verbose=TRUE,...) {
+ function(object,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,emcontrol=em.control(),verbose=TRUE,...) {
fi <- !is.null(fixed)
cr <- !is.null(conrows)
@@ -20,24 +21,21 @@
} else {
if(method=="EM") {
if(constr) {
- warning("EM not applicable for constrained models; optimization method changed to 'donlp'")
- method="donlp"
+ warning("EM not applicable for constrained models; optimization method changed to 'rsolnp'")
+ 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. ")
-
if(method=="EM") {
- if(is.null(em.control$tol)) em.control$tol <- 1e-8
- if(is.null(em.control$crit)) em.control$crit <- "relative"
- if(is.null(em.control$random.start)) em.control$random.start <- FALSE
- object <- em(object,tol=em.control$tol,crit=em.control$crit,random.start=em.control$random.start,verbose=verbose,...)
+ object <- em(object,maxit=emcontrol$maxit,tol=emcontrol$tol,crit=emcontrol$crit,random.start=emcontrol$random.start,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")
Added: pkg/depmixS4/R/em.control.R
===================================================================
--- pkg/depmixS4/R/em.control.R (rev 0)
+++ pkg/depmixS4/R/em.control.R 2010-10-20 11:01:22 UTC (rev 444)
@@ -0,0 +1,3 @@
+em.control <-
+function(maxit=100,tol=1e-8,crit="relative",random.start=TRUE) {
+ return(list(maxit=maxit,tol=tol,crit=crit,random.start=random.start))}
\ No newline at end of file
Property changes on: pkg/depmixS4/R/em.control.R
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/depmixS4/inst/doc/depmixS4.Rnw
===================================================================
--- pkg/depmixS4/inst/doc/depmixS4.Rnw 2010-10-19 15:08:36 UTC (rev 443)
+++ pkg/depmixS4/inst/doc/depmixS4.Rnw 2010-10-20 11:01:22 UTC (rev 444)
@@ -439,14 +439,18 @@
be generated randomly within the EM algorithm by generating random
uniform values for the values of $\gamma_{t}$ in (\ref{eq:Q}) at
iteration 0. Automatic generation of starting values is not available
-for constrained models (see below).
+for constrained models (see below). In the call to \code{fit} below,
+the argument \code{emc=em.control(rand=FALSE)} controls the EM
+algorithm and specifies that random start values should not be
+generated\footnote{As of version 1.0-1, the default is have random
+parameter initialization when using the EM algorithm.}.
\paragraph{Fitting the model and printing results.} The \code{depmix}
function returns an object of class `\code{depmix}' which contains the
model specification, and not a fitted model. Hence, the model needs
to be fitted by calling \code{fit}:
<<>>=
-fm <- fit(mod)
+fm <- fit(mod, emc=em.control(rand=FALSE))
@
The \code{fit} function returns an object of class
@@ -511,7 +515,7 @@
set.seed(1)
mod <- depmix(rt ~ 1, data = speed, nstates = 2, family = gaussian(),
transition = ~ scale(Pacc), instart = runif(2))
-fm <- fit(mod, verbose = FALSE)
+fm <- fit(mod, verbose = FALSE, emc=em.control(rand=FALSE))
@
Note the use of \code{verbose = FALSE} to suppress printing of
@@ -542,7 +546,7 @@
mod <- depmix(list(rt ~ 1,corr ~ 1), data = speed, nstates = 2,
family = list(gaussian(), multinomial("identity")),
transition = ~ scale(Pacc), instart = runif(2))
-fm <- fit(mod,verbose = FALSE)
+fm <- fit(mod, verbose = FALSE, emc=em.control(rand=FALSE))
@
This provides the following fitted model parameters (only the
@@ -616,7 +620,7 @@
mod <- depmix(list(rt ~ 1,corr ~ 1), data = speed, transition = ~ Pacc,
nstates = 2, family = list(gaussian(), multinomial("identity")),
trstart = trst, instart = c(0.99, 0.01))
-fm1 <- fit(mod,verbose = FALSE)
+fm1 <- fit(mod,verbose = FALSE, emc=em.control(rand=FALSE))
@
After this, we use the fitted values from this model to constrain the
@@ -709,7 +713,7 @@
multinomial("identity"), multinomial("identity"),
multinomial("identity")), respstart = runif(24), prior = ~ age,
initdata = balance)
-fm <- fit(mod, verbose = FALSE)
+fm <- fit(mod, verbose = FALSE, emc=em.control(rand=FALSE))
fm
@
Note here that we define a \code{mix} model instead of a \code{depmix}
@@ -978,7 +982,7 @@
<<>>=
mod <- makeDepmix(response = rModels, transition = transition,
prior = inMod, stat = FALSE)
-fm <- fit(mod, verbose = FALSE)
+fm <- fit(mod, verbose = FALSE, emc=em.control(rand=FALSE))
@
Using \code{summary} will print the fitted parameters. Note that the
Modified: pkg/depmixS4/man/depmix.fit.Rd
===================================================================
--- pkg/depmixS4/man/depmix.fit.Rd 2010-10-19 15:08:36 UTC (rev 443)
+++ pkg/depmixS4/man/depmix.fit.Rd 2010-10-20 11:01:22 UTC (rev 444)
@@ -26,14 +26,12 @@
\usage{
\S4method{fit}{depmix}(object, fixed=NULL, equal=NULL, conrows=NULL,
- conrows.upper=0, conrows.lower=0, method=NULL,em.control=list(tol=1e-8,
- crit=c("relative","absolute"),random.start=FALSE),verbose=TRUE,...)
+ conrows.upper=0, conrows.lower=0, method=NULL, emcontrol=em.control(), verbose=TRUE,...)
\S4method{summary}{depmix.fitted}(object,which="all")
\S4method{fit}{mix}(object, fixed=NULL, equal=NULL, conrows=NULL,
- conrows.upper=0, conrows.lower=0, method=NULL,em.control=list(tol=1e-8,
- crit=c("relative","absolute"),random.start=FALSE),verbose=TRUE,...)
+ conrows.upper=0, conrows.lower=0, method=NULL, emcontrol=em.control(), verbose=TRUE,...)
\S4method{summary}{mix.fitted}(object,which="all")
@@ -56,8 +54,8 @@
\item{method}{The optimization method; mostly determined by
constraints.}
- \item{em.control}{Named list with control parameters for the EM
- algorithm (see Control parameters below).}
+ \item{emcontrol}{Named list with control parameters for the EM
+ algorithm (see \code{\link{em.control}}).}
\item{verbose}{Should optimization information be displayed on screen?}
@@ -72,7 +70,7 @@
Models are fitted by the EM algorithm if there are no constraints on
the parameters. Aspects of the EM algorithm can be controlled through
- the \code{control.em} argument (see below).
+ the \code{emcontrol} argument; see details in \code{\link{em.control}}.
Otherwise the general optimizers \code{solnp}, the
default (from package \code{Rsolnp}) or \code{donlp2} (from package
\code{Rdonlp2}) are used which handle general linear (in-)equality
@@ -91,7 +89,7 @@
estimated to be equal. Any integers can be used in this way except 0
and 1, which indicate fixed and free parameters, respectively.
- Using \code{solnp} or \code{donlp2} , a Newton-Raphson scheme is employed
+ Using \code{solnp} (or \code{donlp2}), a Newton-Raphson scheme is employed
to estimate parameters subject to linear constraints by imposing:
bl <= A*x <= bu,
@@ -101,7 +99,7 @@
The \code{conrows} argument is used to specify rows of A directly, and
the conrows.lower and conrows.upper arguments to specify the bounds on
- the constraints. \code{conrows} is a matrix of npar(object) columns
+ the constraints. \code{conrows} must be a matrix of npar(object) columns
and one row for each constraint (a vector in the case of a single
constraint). Examples of these three ways of constraining parameters
are provided below.
@@ -117,35 +115,6 @@
}
-\section{Control parameters for the EM algorithm}{
-
- Aspects of EM algorithm can be controlled by the \code{em.control}
- argument. This named list currently observes the following parameters:
-
- \describe{
- \item{\code{tol}:}{sets the the tolerance level for convergence.}
-
- \item{\code{criterion}:}{sets the convergence criterion to either
- the relative change in the log-likelihood or the absolute change in
- the log-likelihood. The relative likelihood criterion (the
- default) assumes convergence on iteration \eqn{i}{i} when
- \eqn{\frac{\log L(i) - \log L(i-1)}{\log L(i-1)} <
- tol}{\frac{\log L(i) - \log L(i-1)}{\log L(i-1)} <
- tol}. The absolute likelihood criterion assumes convergence
- on iteration \eqn{i}{i} when \eqn{\log L(i) - \log L(i-1) <
- tol}{(logLik(i) - logLik(i-1)) < tol}.}
-
- \item{\code{random.start}:}{is used for a (limited) random
- initialization of the parameters. In particular, if
- \code{random.start=TRUE}, the (posterior) state probabilities are
- randomized at iteration 0 (using a uniform distribution). Random
- initialization is useful when no initial parameters can be given to
- distinguish between the states. It is also useful for repeated
- estimation from different starting values.
- }
- }
-}
-
\value{
\code{fit} returns an object of class
@@ -155,7 +124,7 @@
\describe{
\item{\code{message}:}{Convergence information.}
- \item{\code{conMa}t:}{The constraint matrix A, see Details.}
+ \item{\code{conMat}:}{The constraint matrix A, see Details.}
\item{\code{posterior}:}{The posterior state sequence (computed
with the viterbi algorithm), and the posterior probabilities (delta
@@ -195,14 +164,10 @@
# 2-state model on rt and corr from speed data set
# with Pacc as covariate on the transition matrix
-# starting values for the transition pars (without
-# those EM does not get off the ground)
-set.seed(1)
-tr=runif(6)
-trst=c(tr[1:2],0,tr[3:5],0,tr[6])
mod1 <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,nstates=2,
- family=list(gaussian(),multinomial("identity")),trstart=trst)
+ family=list(gaussian(),multinomial("identity")),ntimes=c(168,134,137))
# fit the model
+set.seed(3)
fmod1 <- fit(mod1)
fmod1 # to see the logLik and optimization information
# to see the parameters
@@ -268,25 +233,18 @@
data(balance)
# four binary items on the balance scale task
-
-instart=c(0.5,0.5)
-set.seed(1)
-respstart=runif(16)
-# note that ntimes argument is used to make this a mixture model
mod4 <- mix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
- family=list(multinomial(),multinomial(),multinomial(),multinomial()),
- respstart=respstart,instart=instart)
+ family=list(multinomial("identity"),multinomial("identity"),multinomial("identity"),multinomial("identity")))
+set.seed(1)
fmod4 <- fit(mod4)
# add age as covariate on class membership by using the prior argument
-instart=c(0.5,0.5,0,0) # we need the initial probs and the coefficients of age
-set.seed(2)
-respstart=runif(16)
mod5 <- mix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
- family=list(multinomial(),multinomial(),multinomial(),multinomial()),
- instart=instart, respstart=respstart, prior=~age, initdata=balance)
+ family=list(multinomial("identity"),multinomial("identity"),multinomial("identity"),multinomial("identity")),
+ prior=~age, initdata=balance)
+set.seed(1)
fmod5 <- fit(mod5)
# check the likelihood ratio; adding age significantly improves the goodness-of-fit
Added: pkg/depmixS4/man/em.control.Rd
===================================================================
--- pkg/depmixS4/man/em.control.Rd (rev 0)
+++ pkg/depmixS4/man/em.control.Rd 2010-10-20 11:01:22 UTC (rev 444)
@@ -0,0 +1,60 @@
+\name{em.control}
+
+\alias{em.control}
+
+\title{Control parameters for the EM algorithm}
+
+\description{Set control parameters for the EM algorithm.}
+
+\usage{
+
+ em.control(maxit = 100, tol = 1e-08, crit = "relative", random.start = TRUE)
+
+}
+
+\arguments{
+
+ \item{maxit}{The maximum number of iterations.}
+
+ \item{tol}{The tolerance level for convergence. See Details.}
+
+ \item{crit}{Sets the convergence criterion to "relative" or "absolute"
+ change of the log-likelihood. See Details.}
+
+ \item{random.start}{This is used for a (limited) random
+ initialization of the parameters. See Details.}
+
+}
+
+\details{
+
+The argument \code{crit} sets the convergence criterion to either the
+relative change in the log-likelihood or the absolute change in the
+log-likelihood. The relative likelihood criterion (the default) assumes
+convergence on iteration \eqn{i}{i} when \eqn{\frac{\log L(i) - \log
+L(i-1)}{\log L(i-1)} < tol}{\frac{\log L(i) - \log L(i-1)}{\log L(i-1)} <
+tol}. The absolute likelihood criterion assumes convergence on iteration
+\eqn{i}{i} when \eqn{\log L(i) - \log L(i-1) < tol}{(logLik(i) -
+logLik(i-1)) < tol}. Use \code{crit="absolute"} to invoke the latter
+convergence criterion. Note that in that case, optimal values of the
+tolerance parameter \code{tol} scale with the value of the log-likelihood.
+
+Argument \code{random.start} This is used for a (limited) random
+initialization of the parameters. In particular, if
+\code{random.start=TRUE}, the (posterior) state probabilities are
+randomized at iteration 0 (using a uniform distribution). Random
+initialization is useful when no initial parameters can be given to
+distinguish between the states. It is also useful for repeated estimation
+from different starting values.
+
+}
+
+\value{
+
+ \code{em.control} returns a list of the control parameters.
+
+}
+
+\author{Ingmar Visser & Maarten Speekenbrink}
+
+\keyword{methods}
\ No newline at end of file
Property changes on: pkg/depmixS4/man/em.control.Rd
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/depmixS4/man/makeDepmix.Rd
===================================================================
--- pkg/depmixS4/man/makeDepmix.Rd 2010-10-19 15:08:36 UTC (rev 443)
+++ pkg/depmixS4/man/makeDepmix.Rd 2010-10-20 11:01:22 UTC (rev 444)
@@ -89,42 +89,46 @@
\examples{
-# below example recreates the model from the depmix help page albeit in a
-# roundabout way
+# below example recreates the same model as on the fit help page in a roundabout way
+# there we had:
+# mod1 <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,nstates=2,
+# family=list(gaussian(),multinomial("identity")),ntimes=c(168,134,137))
data(speed)
rModels <- list(
list(
- GLMresponse(formula=rt~1,data=speed,family=gaussian(),pstart=c(5.52,.202)),
- GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(0.5,0.5))
+ GLMresponse(formula=rt~1,data=speed,family=gaussian()),
+ GLMresponse(formula=corr~1,data=speed,family=multinomial("identity"))
),
list(
- GLMresponse(formula=rt~1,data=speed,family=gaussian(),pstart=c(6.39,.24)),
- GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(.1,.9))
+ GLMresponse(formula=rt~1,data=speed,family=gaussian()),
+ GLMresponse(formula=corr~1,data=speed,family=multinomial("identity"))
)
)
-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]))
+transition[[1]] <- transInit(~Pacc,nstates=2,data=speed)
+transition[[2]] <- transInit(~Pacc,nstates=2,data=speed)
-instart=c(0,1)
-inMod <- transInit(~1,ns=2,ps=instart,data=data.frame(rep(1,3)))
-mod <- makeDepmix(response=rModels,transition=transition,prior=inMod,ntimes=attr(speed,"ntimes"))
+inMod <- transInit(~1,ns=2,data=data.frame(rep(1,3)),family=multinomial("identity"))
+mod <- makeDepmix(response=rModels,transition=transition,prior=inMod,ntimes=c(168,134,137),stationary=FALSE)
+set.seed(3)
fm1 <- fit(mod)
fm1
summary(fm1)
+
# generate data from two different multivariate normal distributions
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)
+set.seed(2)
y1 <- mvrnorm(50,m1,sd1)
y2 <- mvrnorm(50,m2,sd2)
+# this creates data with a single change point
y <- rbind(y1,y2)
# now use makeDepmix to create a depmix model for this bivariate normal timeseries
@@ -143,9 +147,10 @@
mod <- makeDepmix(response=rModels,transition=transition,prior=inMod)
-fm2 <- fit(mod)
+fm2 <- fit(mod,emc=em.control(random=FALSE))
+
# where is the switch point?
-plot(as.ts(posterior(fm2)[,1]))
+plot(as.ts(posterior(fm2)[,2]))
# in below example we add the exgaus distribution as a response model and fit
@@ -284,7 +289,7 @@
mod <- makeDepmix(response=rModels,transition=transition,prior=inMod,ntimes=attr(speed,"ntimes"),stat=FALSE)
-fm3 <- fit(mod)
+fm3 <- fit(mod,emc=em.control(rand=FALSE))
summary(fm3)
}
More information about the depmix-commits
mailing list