[Depmix-commits] r326 - trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 28 17:01:19 CET 2010
Author: ingmarvisser
Date: 2010-01-28 17:01:18 +0100 (Thu, 28 Jan 2010)
New Revision: 326
Modified:
trunk/R/responseMVN.R
Log:
Removed redundant parameters in covariance matrix of multivariate normal models (only the lower triangle of sigma is stored now rather than the whole matrix).
Modified: trunk/R/responseMVN.R
===================================================================
--- trunk/R/responseMVN.R 2010-01-28 11:35:30 UTC (rev 325)
+++ trunk/R/responseMVN.R 2010-01-28 16:01:18 UTC (rev 326)
@@ -10,8 +10,7 @@
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)
object at parameters$coefficients <- fit$coefficients
- if(!is.null(w)) object at parameters$Sigma <- cov.wt(x=fit$residuals,wt=w)$cov else object at parameters$Sigma <- cov(fit$residuals)
- #object <- setpars(object,unlist(pars))
+ if(!is.null(w)) object at parameters$Sigma <- vech(cov.wt(x=fit$residuals,wt=w)$cov) else object at parameters$Sigma <- vech(cov(fit$residuals))
object
}
)
@@ -25,38 +24,38 @@
}
if (missing(mean)) {
mean <- rep(0, length = ncol(x))
- }
- if(missing(invSigma)) {
+ }
+ if(missing(invSigma)) {
if (missing(sigma)) {
sigma <- diag(ncol(x))
- }
- invSigma <- solve(sigma)
- }
- # check consistency
+ }
+ invSigma <- solve(sigma)
+ }
+ # 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(NROW(mean) == NROW(x)) {
- # varying means
-
- # from "mahalanobis":
- x <- as.matrix(x) - mean
- distval <- rowSums((x %*% invSigma) * x)
- #names(retval) <- rownames(x)
- #retval
- } else {
+ }
+ if (NCOL(invSigma) != NCOL(mean)) {
+ stop("mean and sigma have non-conforming size")
+ }
+ if(NROW(mean) == NROW(x)) {
+ # varying means
+
+ # from "mahalanobis":
+ x <- as.matrix(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)
+ }
+ distval <- mahalanobis(x, center = mean, cov = invSigma, inverted=TRUE)
}
if(missing(logdet)) logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
@@ -70,17 +69,17 @@
setMethod("logDens","MVNresponse",
function(object,...) {
- dm_dmvnorm(x=object at y,mean=predict(object),sigma=object at parameters$Sigma,log=TRUE,...)
+ dm_dmvnorm(x=object at y,mean=predict(object),sigma=xpnd(object at parameters$Sigma),log=TRUE,...)
}
)
setMethod("dens","MVNresponse",
function(object,log=FALSE,...) {
- dm_dmvnorm(x=object at y,mean=predict(object),sigma=object at parameters$Sigma,log=log,...)
+ dm_dmvnorm(x=object at y,mean=predict(object),sigma=xpnd(object at parameters$Sigma),log=log,...)
}
-)
-
+)
+
setMethod("predict","MVNresponse",
function(object) {
object at x%*%object at parameters$coefficients
@@ -88,22 +87,28 @@
)
setMethod("simulate",signature(object="MVNresponse"),
- function(object,nsim=1,seed=NULL,times) {
- if(!is.null(seed)) set.seed(seed)
- if(missing(times)) {
- # draw in one go
- mu <- predict(object)
- } else {
- mu <- predict(object)[times,]
- }
- nt <- nrow(mu)
- response <- mvrnorm(nt*nsim,mu=mu,Sigma=object at parameters$Sigma)
- #if(nsim > 1) response <- array(response,dim=c(nt,ncol(response),nsim))
- return(response)
- }
-)
-setMethod("MVNresponse",
- signature(formula="formula"),
+ function(object,nsim=1,seed=NULL,times) {
+ if(!is.null(seed)) set.seed(seed)
+ if(missing(times)) {
+ # draw in one go
+ mu <- predict(object)
+ } else {
+ mu <- predict(object)[times,]
+ }
+ nt <- nrow(mu)
+ if(nrow(object at parameters$coefficients==1)) response <- mvrnorm(nt*nsim,mu=mu[1,],Sigma=xpnd(object at parameters$Sigma))
+ else {
+ response <- matrix(0,nrow(mu),ncol(mu))
+ for(i in 1:nrow(mu)) {
+ response[i,] <- response <- mvrnorm(1,mu=mu[i,],Sigma=xpnd(object at parameters$Sigma))
+ }
+ }
+ return(response)
+ }
+)
+
+setMethod("MVNresponse",
+ signature(formula="formula"),
function(formula,data,pstart=NULL,fixed=NULL,...) {
call <- match.call()
mf <- match.call(expand.dots = FALSE)
@@ -117,23 +122,18 @@
if(!is.matrix(y)) y <- matrix(y,ncol=1)
parameters <- list()
parameters$coefficients <- matrix(0.0,ncol=ncol(y),nrow=ncol(x))
- parameters$Sigma <- diag(ncol(y))
- npar <- length(unlist(parameters))
+ parameters$Sigma <- vech(diag(ncol(y)))
+ npar <- length(unlist(parameters))
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$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
- pstart <- matrix(pstart,ncol(x),byrow=TRUE)
- if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),]
- if(length(unlist(parameters))>length(parameters$coefficients)) {
- tmp <- as.numeric(pstart[(length(parameters$coefficients)+1):length(pstart)])
- if(length(tmp) == ncol(parameters$Sigma)) parameters$Sigma <- diag(tmp) else parameters$Sigma <- matrix(tmp,ncol=ncol(y),nrow=ncol(y))
- }
+ parameters$coefficients <- matrix(pstart[1:(ncol(x)*ncol(y))],ncol(x),byrow=T)
+ parameters$Sigma <- as.numeric(pstart[(length(parameters$coefficients)+1):length(pstart)])
}
mod <- new("MVNresponse",formula=formula,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar)
mod
}
-)
+)
setMethod("show","MVNresponse",
function(object) {
@@ -142,7 +142,7 @@
cat("Coefficients: \n")
print(object at parameters$coefficients)
cat("Sigma: \n")
- print(object at parameters$Sigma)
+ print(xpnd(object at parameters$Sigma))
}
)
@@ -155,10 +155,8 @@
switch(which,
"pars" = {
object at parameters$coefficients <- matrix(values[1:length(object at parameters$coefficients)],ncol(object at x))
- if(length(unlist(object at parameters))>length(object at parameters$coefficients)) {
- st <- length(object at parameters$coefficients)+1
- object at parameters$Sigma <- matrix(as.numeric(values[st:(st+length(object at parameters$Sigma))]),ncol=ncol(object at parameters$Sigma),nrow=nrow(object at parameters$Sigma))
- }
+ st <- length(object at parameters$coefficients)+1
+ object at parameters$Sigma <- as.numeric(values[st:(st+length(object at parameters$Sigma)-1)])
},
"fixed" = {
object at fixed <- as.logical(values)
@@ -168,6 +166,7 @@
return(object)
}
)
+
setMethod("getpars","MVNresponse",
function(object,which="pars",...) {
switch(which,
@@ -182,4 +181,4 @@
)
return(pars)
}
-)
+)
More information about the depmix-commits
mailing list