[Depmix-commits] r668 - pkg/depmixS4/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 11 17:03:34 CEST 2018
Author: maarten
Date: 2018-09-11 17:03:34 +0200 (Tue, 11 Sep 2018)
New Revision: 668
Modified:
pkg/depmixS4/R/EM.R
pkg/depmixS4/R/responseMVN.R
pkg/depmixS4/R/responseNORM.R
Log:
Fixed issue with MVNnorm using unbiased rather than ML estimate of covariance; message provides a little more information when likelihood increases during EM iteration
Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R 2018-08-24 12:30:29 UTC (rev 667)
+++ pkg/depmixS4/R/EM.R 2018-09-11 15:03:34 UTC (rev 668)
@@ -217,7 +217,7 @@
# this should not really happen...
if(j > 0 && (LL.old - LL) >= tol) {
likelihood_decreased <- TRUE
- warning("likelihood decreased on iteration ",j)
+ warning("Log likelihood decreased on iteration ",j, " from ", LL.old, " to ", LL)
break
}
if(j > 0 && ((crit == "absolute" && abs(LL.old - LL) < tol) || (crit == "relative" && abs(LL - LL.old)/abs(LL.old) < tol))) {
@@ -410,7 +410,7 @@
# this should not really happen...
if(j > 0 && (LL.old - LL) > tol) {
likelihood_decreased <- TRUE
- warning("likelihood decreased on iteration ",j)
+ warning("Log likelihood decreased on iteration ",j, " from ", LL.old, " to ", LL)
break
}
if(j > 0 && ((crit == "absolute" && abs(LL.old - LL) < tol) || (crit == "relative" && abs(LL - LL.old)/abs(LL.old) < tol))) {
Modified: pkg/depmixS4/R/responseMVN.R
===================================================================
--- pkg/depmixS4/R/responseMVN.R 2018-08-24 12:30:29 UTC (rev 667)
+++ pkg/depmixS4/R/responseMVN.R 2018-09-11 15:03:34 UTC (rev 668)
@@ -24,67 +24,164 @@
function(object,w) {
if(missing(w)) w <- NULL
pars <- object at parameters
- if(!is.null(w)) fit <- lm.wfit(x=object at x,y=object at y,w=w) else fit <- lm.fit(x=object at x,y=object at y)
+ nas <- any(is.na(object at y))
+ # if(NCOL(object at x) == 1 && all(object at x == 1)) {
+ # # intercept only model, use standard MVN estimation for EM
+ # if(!is.null(w)) {
+ # if(nas) {
+ # y <- object at y[rowSums(is.na(object at y))==0,,drop=FALSE]
+ # w <- w[rowSums(is.na(object at y))==0]
+ # mean <- apply(y,2,weighted.mean,w=w)
+ # cov <- cov.wt(y,wt=w,method = "ML")$cov
+ # } else {
+ # mean <- apply(object at y,2,weighted.mean,w=w)
+ # cov <- cov.wt(object at y,wt=w,method = "ML")$cov
+ # }
+ # } else {
+ # if(nas) {
+ # y <- object at y[rowSums(is.na(object at y))==0,,drop=FALSE]
+ # mean <- colMeans(y)
+ # cov <- cov.wt(y,method="ML")$cov
+ # } else {
+ # mean <- colMeans(object at y)
+ # cov <- cov.wt(object at y,method="ML")$cov
+ # }
+ # }
+ # object at parameters$coefficients <- mean
+ # object at parameters$Sigma <- cov2par(cov)
+ # return(object)
+ # }
+ if(!is.null(w)) {
+ if(nas) {
+ x <- object at x[rowSums(is.na(object at y))==0,,drop=FALSE]
+ y <- object at y[rowSums(is.na(object at y))==0,,drop=FALSE]
+ w <- w[rowSums(is.na(object at y))==0]
+ fit <- lm.wfit(x=x,y=y,w=w)
+ } else {
+ fit <- lm.wfit(x=object at x,y=object at y,w=w) #else fit <- lm.fit(x=object at x,y=object at y)
+ }
+ } else {
+ if(nas) {
+ x <- object at x[rowSums(is.na(object at y))==0,,drop=FALSE]
+ y <- object at y[rowSums(is.na(object at y))==0,,drop=FALSE]
+ fit <- lm.fit(x=x,y=y)
+ } else {
+ fit <- lm.fit(x=object at x,y=object at y)
+ }
+ }
object at parameters$coefficients <- fit$coefficients
- if(!is.null(w)) object at parameters$Sigma <- cov2par(cov.wt(x=fit$residuals,wt=w)$cov) else object at parameters$Sigma <- cov2par(cov(fit$residuals))
- object
+ if(!is.null(w)) {
+ object at parameters$Sigma <- cov2par(cov.wt(x=fit$residuals,wt=w,method="ML")$cov)
+ } else {
+ object at parameters$Sigma <- cov2par(cov.wt(x=fit$residuals,method="ML")$cov)
+ }
+ return(object)
}
)
+
dm_dmvnorm <- function(x,mean,sigma,log=FALSE,logdet,invSigma) {
+ # update 29/08/2018: using new function definition from mvtnorm 1.0-8
# taken from mvtnorm package
+ # using Cholesky decomposition
# allows passing of logdet (sigma) and invsigma to save
- # computation when called repeated times with same sigma
- if (is.vector(x)) {
- x <- matrix(x, ncol = length(x))
+ # computation when called repeated times with same sigma
+ #
+ # mean can be a vector of length ncol(x) or a matrix of dim(x)
+ if(is.vector(x)) {
+ x <- matrix(x, ncol = length(x))
+ }
+ p <- ncol(x)
+ if(!missing(invSigma)) {
+ warning("Use of logdet and invsigma are depreciated")
+ # do the old style computation
+ # as far as I'm aware, depmixS4 doesn't pass invSigma or logdet to this function
+ # so this shouldn't really be necessary
+ if(missing(mean)) {
+ mean <- matrix(0, ncol = p)
}
- if (missing(mean)) {
- mean <- matrix(0, ncol = ncol(x))
- }
if(is.vector(mean)) {
- mean <- matrix(mean, ncol = ncol(x))
+ mean <- matrix(mean, ncol = p)
}
- if(missing(invSigma)) {
- if (missing(sigma)) {
- sigma <- diag(ncol(x))
- }
- invSigma <- solve(sigma)
+ # check consistency
+ if (NCOL(x) != NCOL(invSigma)) {
+ stop("x and sigma have non-conforming size")
}
- # check consistency
- if (NCOL(x) != NCOL(invSigma)) {
- stop("x and sigma have non-conforming size")
- }
- if (NROW(invSigma) != NCOL(invSigma)) {
- stop("sigma must be a square matrix")
- }
- if (NCOL(invSigma) != NCOL(mean)) {
- stop("mean and sigma have non-conforming size")
- }
- if(missing(logdet)) {
- ev <- eigen(sigma, symmetric = TRUE, only.values = TRUE)$values
- if(!all(ev >= 0)) return(rep(NaN,nrow(x))) else logdet <- sum(log(ev))
- }
- if(NROW(mean) == NROW(x)) {
- # varying means
-
- # from "mahalanobis":
- x <- x - mean
- distval <- rowSums((x %*% invSigma) * x)
- #names(retval) <- rownames(x)
- #retval
- } else {
- # constant mean
- if (length(mean) != NROW(invSigma)) {
- stop("mean and sigma have non-conforming size")
- }
- distval <- mahalanobis(x, center = mean, cov = invSigma, inverted=TRUE)
- }
+ if (NROW(invSigma) != NCOL(invSigma)) {
+ stop("sigma must be a square matrix")
+ }
+ if (NCOL(invSigma) != NCOL(mean)) {
+ stop("mean and sigma have non-conforming size")
+ }
+ if(missing(logdet)) {
+ ev <- eigen(sigma, symmetric = TRUE, only.values = TRUE)$values
+ if(!all(ev >= 0)) return(rep(NaN,nrow(x))) else logdet <- sum(log(ev))
+ }
+ if(NROW(mean) == NROW(x)) {
+ # varying means
+
+ # from "mahalanobis":
+ x <- x - mean
+ distval <- rowSums((x %*% invSigma) * x)
+ #names(retval) <- rownames(x)
+ #retval
+ } else {
+ # constant mean
+ if (length(mean) != NROW(invSigma)) {
+ stop("mean and sigma have non-conforming size")
+ }
+ distval <- mahalanobis(x, center = mean, cov = invSigma, inverted=TRUE)
+ }
logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
- if (log) {
- return(logretval)
- } else {
- return(exp(logretval))
+ } else {
+ ## use Cholesky decomposition
+ if(missing(mean)) {
+ mean <- rep(0,p)
}
+ if(missing(sigma)) {
+ sigma <- diag(p)
+ }
+ if(!missing(mean)) {
+ ##if(!is.null(dim(mean)))
+ ##dim(mean) <- NULL
+ if(NCOL(mean) != NCOL(sigma))
+ stop("mean and sigma have non-conforming size")
+ }
+ if (!missing(sigma)) {
+ if(p != ncol(sigma))
+ stop("x and sigma have non-conforming size")
+ if(!isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
+ check.attributes = FALSE))
+ stop("sigma must be a symmetric matrix")
+ }
+ dec <- tryCatch(chol(sigma), error = function(e) e)
+ if (inherits(dec, "error")) {
+ if(all(dim(x) == dim(mean))) {
+ x.is.mu <- rowSums(x != mean) == 0
+ } else {
+ if(length(mean) == p) {
+ x.is.mu <- colSums(t(x) != mean) == 0
+ }
+ }
+ #x.is.mu <- colSums(t(x) != mean) == 0
+ logretval <- rep.int(-Inf, nrow(x))
+ logretval[x.is.mu] <- Inf
+ }
+ else {
+ if(all(dim(x) == dim(mean))) {
+ tmp <- backsolve(dec, t(x - mean), transpose = TRUE)
+ } else {
+ tmp <- backsolve(dec, t(x) - mean, transpose = TRUE)
+ }
+ rss <- colSums(tmp^2)
+ logretval <- -sum(log(diag(dec))) - 0.5 * p * log(2 * pi) - 0.5 * rss
+ }
+ }
+ if(log) {
+ return(logretval)
+ } else {
+ return(exp(logretval))
+ }
}
setMethod("logDens","MVNresponse",
@@ -139,9 +236,18 @@
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
x <- model.matrix(attr(mf, "terms"),mf)
- if(any(is.na(x))) stop("'depmixS4' does not currently handle covariates with missing data.")
y <- model.response(mf)
if(!is.matrix(y)) y <- matrix(y,ncol=1)
+
+ nas <- any(is.na(y))
+ if(nas) {
+ y_missing <- rowSums(is.na(y))
+ x_missing <- rowSums(is.na(x))
+ if(sum(x_missing[y_missing > 0]) > 0) {
+ stop("'depmixS4' does not currently handle covariates with missing data at places where the response is not missing.")
+ }
+ }
+ #if(any(is.na(x)))
parameters <- list()
parameters$coefficients <- matrix(0.0,ncol=ncol(y),nrow=ncol(x))
parameters$Sigma <- cov2par(diag(ncol(y)))
Modified: pkg/depmixS4/R/responseNORM.R
===================================================================
--- pkg/depmixS4/R/responseNORM.R 2018-08-24 12:30:29 UTC (rev 667)
+++ pkg/depmixS4/R/responseNORM.R 2018-09-11 15:03:34 UTC (rev 668)
@@ -19,7 +19,7 @@
if(!is.null(w)) {
pars$sd <- sqrt(sum(w[!nas]*fit$residuals^2/sum(w[!nas])))
} else {
- pars$sd <- sd(fit$residuals)
+ pars$sd <- sqrt(sum(fit$residuals^2)/length(fit$residuals))
}
object <- setpars(object,unlist(pars))
object
More information about the depmix-commits
mailing list