[Depmix-commits] r565 - in pkg/depmixS4: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 28 20:02:50 CEST 2012
Author: maarten
Date: 2012-08-28 20:02:50 +0200 (Tue, 28 Aug 2012)
New Revision: 565
Modified:
pkg/depmixS4/R/EM.R
pkg/depmixS4/R/depmixfit.R
pkg/depmixS4/R/logLik.R
pkg/depmixS4/man/depmix-methods.Rd
pkg/depmixS4/man/depmix.fitted-class.Rd
pkg/depmixS4/man/mix.fitted-class.Rd
Log:
- logLik now has a method="classification" option
- put computation of "posterior" inside em (still inside mix.fit and depmix.fit when using numerical optimization)
- updated help files
Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R 2012-08-28 13:40:57 UTC (rev 564)
+++ pkg/depmixS4/R/EM.R 2012-08-28 18:02:50 UTC (rev 565)
@@ -178,7 +178,13 @@
}
- if(clsf == "hard") class(object) <- "mix.fitted.classLik" else class(object) <- "mix.fitted"
+ if(clsf == "hard") {
+ class(object) <- "mix.fitted.classLik"
+ data.frame(state=viterbi(object)[,1])
+ } else {
+ class(object) <- "mix.fitted"
+ object at posterior <- viterbi(object)
+ }
if(converge) {
if(clsf == "hard") {
@@ -206,7 +212,7 @@
}
# em for hidden markov models
-em.depmix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,classification=c("soft","hard"),...) {
+em.depmix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,classification=c("soft","hard"),na.allow=TRUE,...) {
if(!is(object,"depmix")) stop("object is not of class 'depmix'")
@@ -243,12 +249,14 @@
# initial expectation
if(clsf == "hard") {
- fbo <- list()
+ fbo <- list()
vstate <- viterbi(object)[,1]
fbo$gamma <- as.matrix(model.matrix(~ factor(vstate,levels=1:ns) - 1))
fbo$xi <- array(0,dim=c(sum(ntimes),ns,ns))
fbo$xi[cbind(1:(sum(ntimes)- 1),vstate[-1],vstate[-length(vstate)])] <- 1
- fbo$logLike <- sum(log((apply(object at dens,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+ B <- object at dens
+ if(na.allow) B[is.na(B)] <- 1
+ fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
} else {
fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
}
@@ -264,7 +272,9 @@
fbo$gamma <- as.matrix(model.matrix(~ factor(vstate,levels=1:ns) - 1))
fbo$xi <- array(0,dim=c(sum(ntimes),ns,ns))
fbo$xi[cbind(1:(sum(ntimes)- 1),vstate[-1],vstate[-length(vstate)])] <- 1
- fbo$logLike <- sum(log((apply(object at dens,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+ B <- object at dens
+ if(na.allow) B[is.na(B)] <- 1
+ fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
} else {
fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
}
@@ -318,7 +328,9 @@
fbo$gamma <- as.matrix(model.matrix(~ factor(vstate,levels=1:ns) - 1))
fbo$xi <- array(0,dim=c(sum(ntimes),ns,ns))
fbo$xi[cbind(1:(sum(ntimes)- 1),vstate[-1],vstate[-length(vstate)])] <- 1
- fbo$logLike <- sum(log((apply(object at dens,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+ B <- object at dens
+ if(na.allow) B[is.na(B)] <- 1
+ fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
} else {
# expectation
fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
@@ -342,7 +354,13 @@
}
- if(clsf == "hard") class(object) <- "depmix.fitted.classLik" else class(object) <- "depmix.fitted"
+ if(clsf == "hard") {
+ class(object) <- "depmix.fitted.classLik"
+ object at posterior <- data.frame(state=viterbi(object)[,1])
+ } else {
+ class(object) <- "depmix.fitted"
+ object at posterior <- viterbi(object)
+ }
if(converge) {
if(clsf == "hard") {
Modified: pkg/depmixS4/R/depmixfit.R
===================================================================
--- pkg/depmixS4/R/depmixfit.R 2012-08-28 13:40:57 UTC (rev 564)
+++ pkg/depmixS4/R/depmixfit.R 2012-08-28 18:02:50 UTC (rev 565)
@@ -222,10 +222,12 @@
object at lin.upper <- lin.u
object at lin.lower <- lin.l
+ object at posterior <- viterbi(object)
+
}
- object at posterior <- viterbi(object)
+
return(object)
}
)
Modified: pkg/depmixS4/R/logLik.R
===================================================================
--- pkg/depmixS4/R/logLik.R 2012-08-28 13:40:57 UTC (rev 564)
+++ pkg/depmixS4/R/logLik.R 2012-08-28 18:02:50 UTC (rev 565)
@@ -1,25 +1,58 @@
# depends on getpars and nobs
setMethod("logLik",signature(object="depmix"),
#function(object,method="lystig") {
- function(object,method="fb") { #4/5/2012: set to fb as this is now in C
- if(method=="fb") ll <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=object at ntimes,stationary=object at stationary)$logLike
+ function(object,method=c("fb","lystig","classification"),na.allow=TRUE) {
+ #4/5/2012: set to fb as this is now in C
+ method <- match.arg(method)
+ if(method=="fb") ll <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=object at ntimes,stationary=object at stationary,na.allow=na.allow)$logLike
if(method=="lystig") ll <- lystig(init=object at init,A=object at trDens,B=object at dens,ntimes=object at ntimes,stationary=object at stationary)$logLike
- attr(ll, "df") <- freepars(object)
+ if(method=="classification") {
+ ns <- nstates(object)
+ ntimes <- ntimes(object)
+ vstate <- viterbi(object)[,1]
+ B <- object at dens
+ if(na.allow) B[is.na(B)] <- 1
+ ll <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+ }
+ attr(ll, "df") <- freepars(object)
attr(ll, "nobs") <- nobs(object)
class(ll) <- "logLik"
ll
}
)
+setMethod("logLik",signature(object="depmix.fitted.classLik"),
+ function(object,method=c("classification","fb","lystig"),na.allow=TRUE) {
+ method <- match.arg(method)
+ callNextMethod(object=object,method=method,na.allow=na.allow)
+ }
+)
+
# depends on getpars and nobs
setMethod("logLik",signature(object="mix"),
#function(object,method="lystig") {
- function(object,method="fb") {
+ function(object,method=c("fb","lystig","classification"),na.allow=TRUE) {
+ method <- match.arg(method)
if(method=="fb") ll <- fb(init=object at init,A=matrix(0,1,1),B=object at dens,ntimes=object at ntimes,stationary=TRUE)$logLike
if(method=="lystig") ll <- lystig(init=object at init,A=matrix(0,1,1),B=object at dens,ntimes=object at ntimes,stationary=TRUE)$logLike
+ if(method=="classification") {
+ ntimes <- ntimes(object)
+ gamma <- fb(init=object at init,A=matrix(0,1,1),B=object at dens,ntimes=ntimes,stationary=TRUE)$gamma
+ vstate <- t(apply(gamma,1,ind.max))
+ B <- object at dens
+ if(na.allow) B[is.na(B)] <- 1
+ ll <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes(object)),vstate)]))
+ }
attr(ll, "df") <- freepars(object)
attr(ll, "nobs") <- nobs(object)
class(ll) <- "logLik"
ll
}
)
+
+setMethod("logLik",signature(object="mix.fitted.classLik"),
+ function(object,method=c("classification","fb","lystig"),na.allow=TRUE) {
+ method <- match.arg(method)
+ callNextMethod(object=object,method=method,na.allow=na.allow)
+ }
+)
Modified: pkg/depmixS4/man/depmix-methods.Rd
===================================================================
--- pkg/depmixS4/man/depmix-methods.Rd 2012-08-28 13:40:57 UTC (rev 564)
+++ pkg/depmixS4/man/depmix-methods.Rd 2012-08-28 18:02:50 UTC (rev 565)
@@ -8,6 +8,9 @@
\alias{logLik,depmix-method}
\alias{logLik,mix-method}
+\alias{logLik,depmix.fitted.classLik-method}
+\alias{logLik,mix.fitted.classLik-method}
+
\alias{nobs}
\alias{nobs,depmix-method}
\alias{nobs,mix-method}
@@ -45,9 +48,12 @@
\usage{
- \S4method{logLik}{depmix}(object,method="fb")
- \S4method{logLik}{mix}(object,method="fb")
+ \S4method{logLik}{depmix}(object,method=c("fb","lystig","classification"),na.allow=TRUE)
+ \S4method{logLik}{mix}(object,method=c("fb","lystig","classification"),na.allow=TRUE)
+ \S4method{logLik}{depmix.fitted.classLik}(object,method=c("classification","fb","lystig"),na.allow=TRUE)
+ \S4method{logLik}{mix.fitted.classLik}(object,method=c("classification","fb","lystig"),na.allow=TRUE)
+
\S4method{nobs}{depmix}(object, ...)
\S4method{nobs}{mix}(object, ...)
@@ -74,12 +80,23 @@
the example.}
\item{method}{The log likelihood can be computed by either the forward
- backward algorithm from Rabiner, 1989, or by the method of Lystig and
- Hughes, 2002. The latter is the default as it is faster because in the
- forward backward routine the state and transition smoothed probabilities
- are also computed which are not neccessary for the log likelihood. Those
- smoothed variables, and the forward and backward variables are accessible
- through the \code{\link{forwardbackward}} function.}
+ backward algorithm (Rabiner, 1989), or by the method of Lystig and
+ Hughes, 2002. The former is the default and implemented in a fast
+ C routine. The forward-backward routine also computes the state and transition
+ smoothed probabilities, which are not directly neccessary for the log likelihood.
+ Those smoothed variables, and the forward and backward variables are accessible
+ through the \code{\link{forwardbackward}} function. When method="classification",
+ the classification likelihood is computed, which is the likelihood of the data
+ assuming the state sequence is known and equal to the maximum a posteriori state
+ sequence. The MAP state sequence is available through the \code{\link{viterbi}}
+ function. The classification likelihood is comuted by default when calling the
+ logLik method on an a model fitted by maximising the classification likelihood.}
+
+ \item{na.allow}{Allow missing observations? When set to FALSE,
+ the logLik method will return NA in the presence of missing observations.
+ When set to TRUE, missing values will be ignored when computing the likelihood.
+ When observations are partly missing (when a multivariate observation has missing
+ values on only some of its dimensionis), this may give unexpected results.}
\item{which}{\code{getpars} function: The default "pars" returns a vector
of all parameters of a \code{depmix} object; the alternative value
Modified: pkg/depmixS4/man/depmix.fitted-class.Rd
===================================================================
--- pkg/depmixS4/man/depmix.fitted-class.Rd 2012-08-28 13:40:57 UTC (rev 564)
+++ pkg/depmixS4/man/depmix.fitted-class.Rd 2012-08-28 18:02:50 UTC (rev 565)
@@ -5,8 +5,11 @@
\alias{depmix.fitted}
\alias{depmix.fitted-class}
-\title{Class "depmix.fitted"}
+\alias{depmix.fitted.classLik}
+\alias{depmix.fitted.classLik-class}
+\title{Class "depmix.fitted" (and "depmix.fitted.classLik")}
+
\description{A fitted \code{\link{depmix}} model.}
\section{Slots}{
@@ -58,8 +61,7 @@
\item{\code{lin.upper}}{The upper bounds on the linear constraints.}
- \item{\code{posterior}:}{Posterior (Viterbi) state sequence (not
- implemented currently).}
+ \item{\code{posterior}:}{Posterior (Viterbi) state sequence.}
}
}
@@ -71,9 +73,11 @@
}
\section{Extends}{
- \code{depmix.fitted} extends the \code{depmix} class.
+ \code{depmix.fitted} extends the \code{"\linkS4class{depmix}"} class directly. \code{depmix.fitted.classLik} is
+ similar to \code{depmix.fitted}, the only difference being that the model is fitted
+ by maximising the classification likelihood.
}
-\author{Ingmar Visser}
+\author{Ingmar Visser & Maarten Speekenbrink}
\keyword{classes}
Modified: pkg/depmixS4/man/mix.fitted-class.Rd
===================================================================
--- pkg/depmixS4/man/mix.fitted-class.Rd 2012-08-28 13:40:57 UTC (rev 564)
+++ pkg/depmixS4/man/mix.fitted-class.Rd 2012-08-28 18:02:50 UTC (rev 565)
@@ -3,8 +3,9 @@
\docType{class}
\alias{mix.fitted-class}
+\alias{mix.fitted.classLik-class}
-\title{Class "mix.fitted"}
+\title{Class "mix.fitted" (and "mix.fitted.classLik")}
\description{A fitted \code{\link{mix}} model.}
@@ -48,8 +49,7 @@
\item{\code{lin.upper}}{The upper bounds on the linear constraints.}
- \item{\code{posterior}:}{Posterior (Viterbi) state sequence (not
- implemented currently).}
+ \item{\code{posterior}:}{Posterior (Viterbi) state sequence.}
}
}
@@ -61,10 +61,12 @@
}
\section{Extends}{
- Class \code{"\linkS4class{mix}"}, directly.
+ Class \code{"\linkS4class{mix}"} directly. \code{mix.fitted.classLik} is
+ similar to \code{mix.fitted}, the only difference being that the model is fitted
+ by maximising the classification likelihood.
}
-\author{Ingmar Visser}
+\author{Ingmar Visser & Maarten Speekenbrink}
\examples{
showClass("mix.fitted")
More information about the depmix-commits
mailing list