[Depmix-commits] r564 - pkg/depmixS4/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 28 15:40:57 CEST 2012
Author: maarten
Date: 2012-08-28 15:40:57 +0200 (Tue, 28 Aug 2012)
New Revision: 564
Modified:
pkg/depmixS4/R/EM.R
pkg/depmixS4/R/depmixfit-class.R
Log:
updates voor classification likelihood
Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R 2012-08-13 10:06:42 UTC (rev 563)
+++ pkg/depmixS4/R/EM.R 2012-08-28 13:40:57 UTC (rev 564)
@@ -36,7 +36,7 @@
}
# em for lca and mixture models
-em.mix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,classification=c("soft","hard"),...) {
+em.mix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,classification=c("soft","hard"),na.allow=TRUE,...) {
clsf <- match.arg(classification)
@@ -73,7 +73,9 @@
if(clsf == "hard") {
fbo <- list()
vstate <- apply(gamma,1,which.max)
- 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,matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
}
@@ -85,19 +87,32 @@
# initial expectation
# treat missing values: FIX ME WITH NA.ALLOW OPTION OR SO.
# WHY NOT USE fb TO COMPUTE LL AND GAMMA?
- B <- object at dens
- B <- replace(B,is.na(B) & !is.nan(B),1)
- B <- apply(object at dens,c(1,3),prod)
+ fbo <- fb(init=object at init,A=matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
+
+ #LL <- fbo$logLike
+ #if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
+
+ #B <- object at dens
+ #B <- replace(B,is.na(B) & !is.nan(B),1)
+ #B <- apply(object at dens,c(1,3),prod)
+
if(clsf == "hard") {
- gamma <- t(apply(object at init*B,1,ind.max))
- LL <- sum(log(rowSums(gamma*B)))
- } else {
- gamma <- object at init*B
- LL <- sum(log(rowSums(gamma)))
- gamma <- gamma/rowSums(gamma)
- }
-
+ #gamma <- t(apply(object at init*B,1,ind.max))
+ #LL <- sum(log(rowSums(gamma*B)))
+ fbo$gamma <- t(apply(fbo$gamma,1,ind.max))
+ vstate <- apply(fbo$gamma,1,which.max)
+ 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)]))
+
+ #LL <- sum(log(rowSums(gamma*B)))
+ } #else {
+ #gamma <- object at init*B
+ #LL <- sum(log(rowSums(gamma)))
+ #gamma <- gamma/rowSums(gamma)
+ #}
+ LL <- fbo$logLike
if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
}
@@ -110,13 +125,13 @@
# should become object at prior <- fit(object at prior)
- object at prior@y <- gamma[bt,,drop=FALSE]
+ object at prior@y <- fbo$gamma[bt,,drop=FALSE]
object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
object at init <- dens(object at prior)
for(i in 1:ns) {
for(k in 1:nresp(object)) {
- object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=gamma[,i])
+ object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=fbo$gamma[,i])
# update dens slot of the model
object at dens[,k,i] <- dens(object at response[[i]][[k]])
}
@@ -125,19 +140,24 @@
# expectation
# treat missing values: FIX ME WITH NA.ALLOW OPTION OR SO.
# WHY NOT USE fb TO COMPUTE LL AND GAMMA?
- B <- object at dens
- B <- replace(B,is.na(B) & !is.nan(B),1)
- B <- apply(object at dens,c(1,3),prod)
+ #B <- object at dens
+ #B <- replace(B,is.na(B) & !is.nan(B),1)
+ #B <- apply(object at dens,c(1,3),prod)
+ fbo <- fb(init=object at init,A=matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
if(clsf == "hard") {
- gamma <- t(apply(object at init*B,1,ind.max))
- LL <- sum(log(rowSums(gamma*B)))
- } else {
- fbo <- fb(init=object at init,matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
- gamma <- fbo$gamma
- LL <- fbo$logLike
- }
-
+ fbo$gamma <- t(apply(fbo$gamma,1,ind.max))
+ vstate <- apply(fbo$gamma,1,which.max)
+ 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)]))
+ #LL <- sum(log(rowSums(gamma*B)))
+ } #else {
+ # fbo <- fb(init=object at init,matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
+ # gamma <- fbo$gamma
+ # LL <- fbo$logLike
+ #}
+ LL <- fbo$logLike
# print stuff
if(verbose&((j%%5)==0)) {
cat("iteration",j,"logLik:",LL,"\n")
@@ -158,7 +178,7 @@
}
- class(object) <- "mix.fitted"
+ if(clsf == "hard") class(object) <- "mix.fitted.classLik" else class(object) <- "mix.fitted"
if(converge) {
if(clsf == "hard") {
@@ -267,7 +287,6 @@
trm <- matrix(0,ns,ns)
for(i in 1:ns) {
if(!object at stationary) {
-
object at transition[[i]]@y <- fbo$xi[,,i]/fbo$gamma[,i]
object at transition[[i]] <- fit(object at transition[[i]],w=as.matrix(fbo$gamma[,i]),ntimes=ntimes(object)) # check this
} else {
@@ -323,7 +342,7 @@
}
- class(object) <- "depmix.fitted"
+ if(clsf == "hard") class(object) <- "depmix.fitted.classLik" else class(object) <- "depmix.fitted"
if(converge) {
if(clsf == "hard") {
Modified: pkg/depmixS4/R/depmixfit-class.R
===================================================================
--- pkg/depmixS4/R/depmixfit-class.R 2012-08-13 10:06:42 UTC (rev 563)
+++ pkg/depmixS4/R/depmixfit-class.R 2012-08-28 13:40:57 UTC (rev 564)
@@ -20,6 +20,17 @@
contains="mix"
)
+setClass("mix.fitted.classLik",
+ representation(
+ message="character", # convergence information
+ conMat="matrix", # constraint matrix on the parameters for general linear constraints
+ lin.upper="numeric", # upper bounds for linear constraint
+ lin.lower="numeric", # lower bounds for linear constraints
+ posterior="data.frame" # posterior probabilities for the states
+ ),
+ contains="mix"
+)
+
# accessor functions
setMethod("posterior","mix.fitted",
@@ -56,6 +67,15 @@
contains="depmix"
)
+setClass("depmix.fitted.classLik",
+ representation(message="character", # convergence information
+ conMat="matrix", # constraint matrix on the parameters for general linear constraints
+ lin.upper="numeric", # upper bounds for linear constraints
+ lin.lower="numeric", # lower bounds for linear constraints
+ posterior="data.frame" # posterior probabilities for the states
+ ),
+ contains="depmix"
+)
# accessor functions
setMethod("posterior","depmix.fitted",
More information about the depmix-commits
mailing list