[Depmix-commits] r293 - papers/jss
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 13 14:36:57 CEST 2009
Author: ingmarvisser
Date: 2009-07-13 14:36:51 +0200 (Mon, 13 Jul 2009)
New Revision: 293
Modified:
papers/jss/article.tex
Log:
Added all the code examples, albeit with minimal text.
Modified: papers/jss/article.tex
===================================================================
--- papers/jss/article.tex 2009-07-10 09:40:28 UTC (rev 292)
+++ papers/jss/article.tex 2009-07-13 12:36:51 UTC (rev 293)
@@ -15,7 +15,8 @@
\Shorttitle{depmixS4 Hidden Markov Models} %% a short title (if necessary)
%% an abstract and keywords
-\Abstract{
+ \Abstract{
+
\pkg{depmixS4} implements a general framework for definining and
fitting hidden Markov mixture model in the R programming language
\citep{R2009}. This includes standard Markov models,
@@ -24,13 +25,16 @@
multivariate data with multinomial and/or gaussian distributions.
Parameters can be estimated subject to general linear constraints.
Parameter estimation is done through an EM algorithm or by a
- direct optimization approach using the
- \pkg{Rdonlp2} optimization routine when contraints are imposed on
- the parameters. A number of illustrative examples are included.
+ direct optimization approach using the \pkg{Rdonlp2} optimization
+ routine when contraints are imposed on the parameters. A number
+ of illustrative examples are included.
}
+
\Keywords{hidden Markov model, dependent mixture model, mixture model, \proglang{R}}
+
\Plainkeywords{hidden Markov model, dependent mixture model, mixture model,R} %% without formatting
+
%% at least one keyword must be supplied
%% publication information
@@ -86,15 +90,16 @@
those latter areas of application, latent Markov models are usually
referred to as hidden Markov models.
-The \pkg{depmixS4} package was motivated by the fact that Markov models
-are used commonly in the social sciences, but no comprehensive package
-was available for fitting such models. Common programs for Markovian
-models include Panmark \citep{Pol1996}, and for latent class models
-Latent Gold \citep{Vermunt2003}. Those programs are lacking a number
-of important features, besides not being freely available. There are
-currently some packages in R that handle hidden Markov models but they
-lack a number of features that we needed in our research. In
-particular, \pkg{depmixS4} was designed to meet the following goals:
+The \pkg{depmixS4} package was motivated by the fact that Markov
+models are used commonly in the social sciences, but no comprehensive
+package was available for fitting such models. Common programs for
+Markovian models include Panmark \citep{Pol1996}, and for latent class
+models Latent Gold \citep{Vermunt2003}. Those programs are lacking a
+number of important features, besides not being freely available.
+There are currently some packages in R that handle hidden Markov
+models but they lack a number of features that we needed in our
+research. In particular, \pkg{depmixS4} was designed to meet the
+following goals:
\begin{enumerate}
\item to be able to handle parameter estimates subject to general
@@ -114,16 +119,16 @@
\end{enumerate}
-Although \pkg{depmixS4} is designed to deal with
-longitudinal or time series data, for say $T>100$, it can also handle
-the limit case with $T=1$. In those cases, there are no time
-dependencies between observed data, and the model reduces to a finite
-mixture model, or a latent class model. Although there are other
-specialized packages to deal with mixture data, one specific feature
-that we needed ourselves which is to the best of our knowledge not
-available in other packages is the possibility to include covariates
-on the prior probabilities of class membership. In the next section,
-an outline is provided of the model and the likelihood equations.
+Although \pkg{depmixS4} is designed to deal with longitudinal or time
+series data, for say $T>100$, it can also handle the limit case with
+$T=1$. In those cases, there are no time dependencies between
+observed data, and the model reduces to a finite mixture model, or a
+latent class model. Although there are other specialized packages to
+deal with mixture data, one specific feature that we needed ourselves
+which is to the best of our knowledge not available in other packages
+is the possibility to include covariates on the prior probabilities of
+class membership. In the next section, an outline is provided of the
+model and the likelihood equations.
\subsection*{Acknowledgements}
@@ -142,19 +147,19 @@
\section{The dependent mixture model}
The data considered here, has the general form $O_{1}^{1}, \ldots,
-O_{1}^{m}$, $O_{2}^{1}, \ldots, O_{2}^{m}$, \ldots, $O_{T}^{1}, \ldots,
-O_{T}^{m}$ for an $m$-variate time series of length $T$. As an
-example, consider a time series of responses generated by a single
+O_{1}^{m}$, $O_{2}^{1}, \ldots, O_{2}^{m}$, \ldots, $O_{T}^{1},
+\ldots, O_{T}^{m}$ for an $m$-variate time series of length $T$. As
+an example, consider a time series of responses generated by a single
subject in a reaction time experiment. The data consists of three
variables, reaction time, accuracy and a covariate which is a pay-off
factor which determines the reward for speed and accuracy. These
-variables are measured on 168, 134 and 137 occasions respectively
-(in Figure~\ref{fig:speed} the first part of this series is plotted).
+variables are measured on 168, 134 and 137 occasions respectively (in
+Figure~\ref{fig:speed} the first part of this series is plotted).
\begin{figure}[htbp]
\begin{center}
- \scalebox{1}{\includegraphics*{graphs/speed1.pdf}}
- \caption{Reaction times, accuracy and pay-off values for
+ \scalebox{1.2}{\includegraphics*{graphs/speed1.pdf}}
+ \caption{Response times, accuracy and pay-off values for
the first series of responses in dataset \texttt{speed}.}
\label{fig:speed}
\end{center}
@@ -169,10 +174,10 @@
likelihood and estimation procedure for the hidden Markov model is
described, given data of the above form. These models are called
dependent mixture models because one of the authors (Ingmar Visser)
-thought it was time for a new name for these models\footnote{Only later
-did I find out that \citet{Leroux1992} already coined the term dependent
-mixture models in an application with hidden Markov mixtures of Poisson
-count data.}
+thought it was time for a new name for these models\footnote{Only
+later did I find out that \citet{Leroux1992} already coined the term
+dependent mixture models in an application with hidden Markov mixtures
+of Poisson count data.}
The dependent mixture model is defined by the following elements:
\begin{enumerate}
@@ -190,12 +195,26 @@
\item a vector $\pmb{\pi}$ of latent state initial probabilities
$\pi_{i}$
\end{enumerate}
-When transitions are added to the latent class model, it is more
-appropriate to refer to the classes as states. The word class is
-rather more associated with a stable trait-like attribute whereas a
+When transitions are added to the latent class model, it is more
+appropriate to refer to the classes as states. The word class is
+rather more associated with a stable trait-like attribute whereas a
state can change over time.
-Need model equations: maybe from Fruhwirth-Schnatter?
+The dependent mixture model is then given by the following equations:
+\begin{align}
+ S_{t} = A S_{t-1}, t=2, \ldots, T \\
+ O_{t} = f(O_{t}|S_{t}),
+\end{align}
+where $S_{t}$ is a sequence of hidden states, $A$ is a transition
+matrix, $O_{t}$ is an (possibly multivariate) observation and $f$ is a
+density function for $O_{t}$ conditional on the hidden state $S_{t}$.
+In the example data above, $f$ could be a gaussian distribution
+function for the response time variable, and a Bernouilli for the
+accuracy data. In the models we are considering here, both the
+transition probabilities $A$ and the initial state probabilities $\pi$
+may depend on covariates as well as the response distributions $f$.
+See for example \citet{Fruhwirth-Schnatter2006} for an overview of
+hidden Markov models with extensions.
\subsection{Likelihood}
@@ -262,6 +281,9 @@
\cite{Tamura2007,Spellucci2002}, because it handles general linear
(in)equality constraints, and optionally also non-linear constraints.
+
+
+
Need some more on EM and how/why it is justified to do separate weighted
fits of the response models and transition and prior models.
@@ -281,37 +303,44 @@
\subsection{Example data: speed}
Throughout this manual a data set called \code{speed} is used. It
-consists of three time series with three variables: reaction time,
+consists of three time series with three variables: response time,
accuracy, and a covariate Pacc which defines the relative pay-off for
-speeded and accurate responding. The participant in this experiment
-switches between fast responding at chance level and relatively
-slower responding at a high level of accuracy.
+speeded and accurate responding. The participant in this experiment
+switches between fast responding at chance level and relatively slower
+responding at a high level of accuracy. Interesting hypotheses to
+test are: is the switching regime symmetric? Is there evidence for
+two states or does one state suffice? Is the guessing state actually
+a guessing state, i.e., is the probability of a correct response at
+chance level (0.5)?
-Interesting hypotheses to test are: is the switching regime symmetric?
-Is there evidence for two states or does one state suffice? Is the
-guessing state actually a guessing state, i.e., is the probability
-of a correct response at chance level (0.5)?
-
-
\subsection{A simple example}
A dependent mixture model is defined by the number of states, and by
the response distribution functions, and can be created with the
-\code{depmix}-function as follows (see the help files for other
-options):
+\code{depmix}-function as follows:
\begin{verbatim}
set.seed(1)
mod <- depmix(rt~1, data=speed, nstates=2, trstart=runif(4))
+\end{verbatim}
+The \code{depmix} function returns an object of class \code{depmix}
+which contains the model specification (and not a fitted model!).
+Note also that start values for the transition parameters are provided
+in this call using the \code{trstart} argument. The package does not
+provide automatic starting values.
+
+The so-defined models needs to be \code{fit}ted with the following:
+\begin{verbatim}
fm <- fit(mod)
\end{verbatim}
-The function returns an object of class \code{depmix.fitted} which
-extends the \code{depmix} class, adding convergence information (and
-information about constraints if these were applied, see section on Constraints below). The
-\code{print} method provides summary information on convergence, the log likelihood
-and the AIC and BIC values. These statistics may be extracted using \code{logLik},
-\code{AIC} and \code{BIC}, respectively.
+The \code{fit}-function returns an object of class
+\code{depmix.fitted} which extends the \code{depmix} class, adding
+convergence information (and information about constraints if these
+were applied, see below). The \code{print} method provides summary
+information on convergence, the log likelihood and the AIC and BIC
+values. These statistics may be extracted using \code{logLik},
+\code{AIC} and \code{BIC}, respectively.
\begin{verbatim}
> fm
@@ -321,8 +350,9 @@
BIC: 211.275
\end{verbatim}
-The \code{summary} method of fitted models provides the parameter estimates, first for
-the prior probabilities model, second for the transition model, and third for the response models.
+The \code{summary} method of \code{fit}ted models provides the parameter
+estimates, first for the prior probabilities model, second for the
+transition model, and third for the response models.
\begin{verbatim}
> summary(fm)
@@ -366,70 +396,188 @@
\end{verbatim}
-\subsection{Adding covariates on transition parameters}
+\subsection{Transition parameters}
The transition matrix is parametrized as a list of multinomial
logistic models. The initial state probabilities are similarly
-parametrized as a multinomial logistic model. Both models use a baseline
-category parametrization, meaning that the parameter for the base
-category is fixed at zero. The default base category is the first
-state. Hence, for example, for a 3-state model, the initial state
-probability model would have three parameters of which the first is
-fixed at zero and the other two are freely estimated.
+parametrized as a multinomial logistic model. Both models use a
+baseline category parametrization, meaning that the parameter for the
+base category is fixed at zero (see \citet[see][p.\ 267
+ff.]{Agresti2002} for multinomial logistic models and various
+parameterizations). The default base category is the first state.
+Hence, for example, for a 3-state model, the initial state probability
+model would have three parameters of which the first is fixed at zero
+and the other two are freely estimated.
-See \citet[see][p.\ 267 ff.]{Agresti2002} for multinomial logistic models and various
-parameterizations.
-
-Covariates can be specified using a one-sided formula as in the
-following example:
-
+Covariates on the transition probabilities can be specified using a
+one-sided formula as in the following example:
\begin{verbatim}
set.seed(1)
-mod <- depmix(list(rt~1,corr~1), data=speed, nstates=2, family=list(gaussian(), multinomial()),
-transition=~Pacc, instart=runif(2))
+mod <- depmix(rt~1, data=speed, nstates=2,
+family=gaussian(), transition=~scale(Pacc),
+instart=runif(2))
fm <- fit(mod)
\end{verbatim}
+\begin{verbatim}
+ Initial state probabilties model
+ Model of type multinomial, formula: ~1
+ Coefficients:
+ [,1] [,2]
+ [1,] 0 10.71779
+ Probalities at zero values of the covariates.
+ 2.214681e-05 0.9999779
+ Transition model for state (component) 1
+ Model of type multinomial, formula: ~scale(Pacc)
+ Coefficients:
+ [,1] [,2]
+ [1,] 0 -0.9215182
+ [2,] 0 1.8649734
+ Probalities at zero values of the covariates.
+ 0.7153513 0.2846487
+ Transition model for state (component) 2
+ Model of type multinomial, formula: ~scale(Pacc)
+ Coefficients:
+ [,1] [,2]
+ [1,] 0 2.471442
+ [2,] 0 3.570856
+ Probalities at zero values of the covariates.
+ 0.07788458 0.9221154
+ Response model(s) for state 1
-\subsection{Add covariates to prior model}
+ Response model for response 1
+ Model of type gaussian, formula: rt ~ 1
+ Coefficients:
+ [1] 5.512179
+ sd 0.1921098
-Add example of balance data.
+ Response model(s) for state 2
-Note that this can be done for the initial state probabilities by
-specifying \code{prior=~X1}, where X1 is the desired covariate. The result
-of this is that the transition probabilities are now dependent on the
-covariate Pacc (which is an experimenter controlled variable to induce
-switching between guessing and accurate responding).
+ Response model for response 1
+ Model of type gaussian, formula: rt ~ 1
+ Coefficients:
+ [1] 6.3885
+ sd 0.2402693
+\end{verbatim}
+
+The summary provides all parameters of the model, also the
+(redundant) zeroes for the base-line category in the multinomial model.
+The summary also prints the transition probabilities at the zero value
+of the covariate. Note that scaling of the covariate is useful in this
+regard as it makes interpretation of these intercept probabilities
+easier.
+
+\subsection{Multivariate data}
+
+Multivariate data can be modelled by providing a list of formulae as
+well as a list of family objects for the distributions of the various
+responses. In above examples we have only used the response times
+which were modelled with the gaussian distribution. The accuracy data
+are in the \code{speed} data are modelled with a multinomial by
+specifying the following:
\begin{verbatim}
-data(balance)
set.seed(1)
-mod <- mix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
- family=list(multinomial(),multinomial(),multinomial(),multinomial()),
- respstart=runif(16))
+mod <- depmix(list(rt~1,corr~1), data=speed, nstates=2,
+ family=list(gaussian(),multinomial()), transition=~scale(Pacc),
+ instart=runif(2))
fm <- fit(mod)
\end{verbatim}
+which provides the following fitted model parameters (only the
+response parameters are given here):
+\begin{verbatim}
+ > summary(fm)
+ \ldots
+ Response model(s) for state 1
+ Response model for response 1
+ Model of type gaussian, formula: rt ~ 1
+ Coefficients:
+ [1] 5.52169
+ sd 0.2028857
+ Response model for response 2
+ Model of type multinomial, formula: corr ~ 1
+ Coefficients:
+ [,1] [,2]
+ [1,] 0 0.1030554
+ Probalities at zero values of the covariates.
+ 0.4742589 0.5257411
+ Response model(s) for state 2
+ Response model for response 1
+ Model of type gaussian, formula: rt ~ 1
+ Coefficients:
+ [1] 6.39369
+ sd 0.2373650
-\subsection{Fixing and constraining parameters}
+ Response model for response 2
+ Model of type multinomial, formula: corr ~ 1
+ Coefficients:
+ [,1] [,2]
+ [1,] 0 2.245514
+ Probalities at zero values of the covariates.
+ 0.09573715 0.9042629
+\end{verbatim}
+As can be seen, state 1 has fast response times around chance level,
+whereas state 2 corresponds with slower responding at higher accuracy
+levels.
-Using package \pkg{Rdonlp2}, parameters may be fitted subject to
-general linear (in-)equality constraints. In particular, the following
-constraints may be applied:
+\subsection{Adding covariates on the prior probabilities}
+The \pkg{depmixS4} contains another data set which is used to
+illustrate the use of covariates on the prior probabilities of models,
+in this case a latent class model. The \code{balance} data consists
+of 4 binary items (correct-incorrect) on a balance scale task
+\citet{Siegler198?}. The data form a subset of the data published in
+\citet{Jansen199?}.
-Constraining and fixing parameters is done using the \code{conpat}
-argument to the \code{depmix.fit}-function, which specifies for each
-parameter in the model whether it's fixed (0) or free (1 or higher).
-Equality constraints can be imposed by having two parameters have the
-same number in the \code{conpat} vector. When only fixed values are
+Similarly to the transition matrix, covariates on the prior
+probabilities of the latent states (or classes in this case), are
+defined by using a one-sided formula:
+\begin{verbatim}
+balance$age <- balance$age-5
+set.seed(1)
+mod <- mix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
+ family=list(multinomial(), multinomial(), multinomial(),
+ multinomial()), respstart=c(rep(c(0.9,0.1),4),rep(c(0.1,0.9),4)),
+ prior=~age, initdata=balance)
+fm <- fit(mod)
+\end{verbatim}
+Note here that we define a \code{mix} model instead of a \code{depmix}
+models as these data form independent observations.
+
+The summary of the \code{fit}ted model gives (only the prior model is
+shown here):
+\begin{verbatim}
+> summary(fm)
+Mixture probabilities model
+Model of type multinomial, formula: ~age
+ [,1] [,2]
+[1,] 0 -2.5182573
+[2,] 0 0.5512996
+Probalities at zero values of the covariates.
+0.9254119 0.07458815
+\end{verbatim}
+Hence at young ages, children have a high probability of incorrect
+responding in class~1, whereas the prior probability for class~2
+increases with age.
+
+
+\subsection{Fixing and constraining parameters}
+
+Using package \pkg{Rdonlp2}, parameters may be fitted subject to
+general linear (in-)equality constraints. Constraining and fixing
+parameters is done using the \code{conpat} argument to the
+\code{depmix.fit}-function, which specifies for each parameter in the
+model whether it's fixed (0) or free (1 or higher). Equality
+constraints can be imposed by having two parameters have the same
+number in the \code{conpat} vector. When only fixed values are
required the \code{fixed} argument can be used instead of
\code{conpat}, with zeroes for fixed parameters and other values (ones
e.g.) for non-fixed parameters. Fitting the models subject to these
@@ -442,38 +590,185 @@
Parameters are numbered in the following order:
\begin{enumerate}
\item the prior model parameters
-
\item the parameters for the transition models
-
\item the response model parameters per state (and subsequently
per response in the case of multivariate time series)
-
\end{enumerate}
To see the ordering of parameters use the following:
\begin{verbatim}
-mod <- setpars(mod, value=1:npar(mod))mod
+mod <- setpars(mod, value=1:npar(mod))
+mod
\end{verbatim}
To see which parameters are fixed (by default only baseline parameters
in the multinomial logistic models for the transition models and the
-initial state probabilities model:
+initial state probabilities model):
\begin{verbatim}
-mod <- setpars(mod,
-getpars(mod,which="fixed"))mod
+mod <- setpars(mod, getpars(mod,which="fixed"))
+mod
\end{verbatim}
+When fitting constraints it is useful to have good starting values
+for the parameters and hence we first fit the following model:
+\begin{verbatim}
+trst=c(0.9,0.1,0,0,0.1,0.9,0,0)
+mod <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,
+ nstates=2,family=list(gaussian(),multinomial()),
+ trstart=trst,inst=c(.999,0.001))
+fm <- fit(mod)
+\end{verbatim}
+After this, we use the fitted values from this model to constrain the
+regression coefficients on the transition matrix:
+\begin{verbatim}
+# start with fixed and free parameters
+conpat <- c(0,1,rep(c(0,1),4),1,1,0,1,1,1,0,1)
+# constrain the beta's on the transition parameters to be equal
+conpat[6] <- conpat[10] <- 2
+fm <- fit(fm,equal=conpat)
+\end{verbatim}
+Using \code{summary} on the fitted model shows that the regression
+coefficients are now estimated at the same value of 12.66.
+Note that these arguments provide the possibility for arbitrary
+constraints, also between, e.g., a multinomial regression coefficient
+for the transition matrix and the mean of a gaussian response model.
+Whether such constraints make sense is hence the responsibility of
+the user.
\section{Extending \pkg{depmixS4}}
-Introduce the exgaus example.
+The \pkg{depmixS4} package was built with the aim of having the
+possibility of adding new response distributions (and possibly also
+other models for the prior and transition probabilities models).
+Therefore, the EM algorithm simply calls \code{fit} functions from the
+separate response models without knowing what they are. As a
+consequence, adding user-specified response models is straightforward.
+User-defined distributions should extend the \code{response}-class
+and have the following slots:
+\begin{enumerate}
+ \item y: the response variable
+ \item x: the design matrix, possibly only an intercept
+ \item paramaters: a named list with the coefficients and possibly
+ other parameters, e.g., the standard deviation in the gaussian
+ response model
+ \item fixed: a vector of logicals indicating whether parameters are
+ fixed
+ \item npar: numerical indicating the number of parameters of the model
+\end{enumerate}
-Explain the use of makeDepmix for having full control of every aspect of the model.
+In the \code{speed} data example, it may be more appropriate to model
+the response times as an exgaus distribution rather than using the
+gaussian. We can do so as follows, first we define our own
+exgaus-class extending the response-class and set it to be a generic:
+\begin{verbatim}
+ setClass("exgaus", contains="response")
+ setGeneric("exgaus", function(y, pstart = NULL, fixed = NULL, ...)
+ standardGeneric("exgaus"))
+\end{verbatim}
+The so-defined class now needs a number of methods:
+\begin{enumerate}
+ \item constructor: function to create instances of the class
+ with starting values
+ \item show method: to print the model to the terminal
+ \item dens: the function that computes the density of the responses
+ \item getpars and setpars: to get and set parameters
+ \item predict: to generate predict'ed values
+ \item fit: function to fit the model using posterior weights (used
+ by the EM algorithm)
+\end{enumerate}
+Only the constructor function and the fit function are provided here
+and the complete code is provided in the helpfile for \code{makeDepmix}.
+The \code{exgaus} example uses the \pkg{gamlss.distr} package for
+computing the \code{dens}ity and for fitting the parameters.
+The constructor function is defined as:
+\begin{verbatim}
+setMethod("exgaus",
+ signature(y="ANY"),
+ function(y,pstart=NULL,fixed=NULL, ...) {
+ y <- matrix(y,length(y))
+ x <- matrix(1)
+ parameters <- list()
+ npar <- 3
+ if(is.null(fixed)) fixed <- as.logical(rep(0,npar))
+ if(!is.null(pstart)) {
+ if(length(pstart)!=npar) stop("length of 'pstart' must be ",npar)
+ parameters$mu <- pstart[1]
+ parameters$sigma <- log(pstart[2])
+ parameters$nu <- log(pstart[3])
+ }
+ mod <- new("exgaus",parameters=parameters,fixed=fixed,x=x,y=y,npar=npar)
+ mod
+ }
+)
+\end{verbatim}
+The \code{fit} function is defined as follows:
+\begin{verbatim}
+setMethod("fit","exgaus",
+ function(object,w) {
+ if(missing(w)) w <- NULL
+ y <- object at y
+ fit <- gamlss(y~1,weights=w,family=exGAUS(),
+ control=gamlss.control(n.cyc=100,trace=FALSE),
+ mu.start=object at parameters$mu,
+ sigma.start=exp(object at parameters$sigma),
+ nu.start=exp(object at parameters$nu))
+ pars <- c(fit$mu.coefficients,fit$sigma.coefficients,fit$nu.coefficients)
+ object <- setpars(object,pars)
+ object
+ }
+)
+The \code{fit} function defines a trivial \code{gamlss} model with
+only an intercept to be estimated and then sets the fitted parameters
+back into their respective slots in the `exgaus' object. See the help
+for \pkg{gamlss.distr} for interpretation of these parameters.
+
+After defining all the necessary methods for the new response model,
+we can now define the dependent mixture model using this reponse.
+The function \code{makeDepmix} is added to \pkg{depmixS4} to have
+full control over model specification, and we need it here.
+
+We first create all the response models that we want as a double list:
+\begin{verbatim}
+rModels <- list(
+ list(
+ exgaus(rt,pstart=c(5,.1,.1)),
+ GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(0.5,0.5))
+ ),
+ list(
+ exgaus(rt,pstart=c(6,.1,.1)),
+ GLMresponse(formula=corr~1,data=speed,family=multinomial(),pstart=c(.1,.9))
+ )
+)
+\end{verbatim}
+Next, we define the transition and prior probability models using the
+transInit function (which produces a transInit model, which also extends
+the response class):
+\begin{verbatim}
+trstart=c(0.9,0.1,0.1,0.9)
+transition <- list()
+transition[[1]] <- transInit(~Pacc,nstates=2,data=speed,pstart=c(0.9,0.1,0,0))
+transition[[2]] <- transInit(~Pacc,nstates=2,data=speed,pstart=c(0.1,0.9,0,0))
+inMod <- transInit(~1,ns=2,pstart=c(0.1,0.9),data=data.frame(1))
+\end{verbatim}
+Finally, we put everything together using \code{makeDepmix} and fit
+the model:
+\begin{verbatim}
+mod <- makeDepmix(response=rModels,transition=transition,
+prior=inMod,stat=FALSE)
+fm <- fit(mod)
+\end{verbatim}
+
+Using \code{summary} will print the fitted parameters. Note that the
+use of \code{makeDepmix} allows the possibility of, say, fitting a
+gaussian in one state and an exgaus distribution in another state.
+
+
+
\section{Conclusion \& future work}
More information about the depmix-commits
mailing list