[Depmix-commits] r366 - in trunk: . R man

[noreply at r-forge.r-project.org]

Author: maarten
Date: 2010-02-24 15:49:50 +0100 (Wed, 24 Feb 2010)
New Revision: 366

Modified:
   trunk/DESCRIPTION
   trunk/R/responseMVN.R
   trunk/man/depmix-internal.Rd
Log:
- added cov2par and par2cov functions for MVNresponse
- removed dependence on MCMpack

Modified: trunk/DESCRIPTION
===================================================================
--- trunk/DESCRIPTION	2010-02-24 10:58:56 UTC (rev 365)
+++ trunk/DESCRIPTION	2010-02-24 14:49:50 UTC (rev 366)
@@ -4,7 +4,7 @@
 Title: Dependent Mixture Models
 Author: Ingmar Visser <i.visser at uva.nl>, Maarten Speekenbrink <m.speekenbrink at ucl.ac.uk>
 Maintainer: Ingmar Visser <i.visser at uva.nl>
-Depends: R (>= 2.9.1), stats, nnet, methods, MASS, MCMCpack
+Depends: R (>= 2.9.1), stats, nnet, methods, MASS
 Suggests: Rdonlp2, gamlss, gamlss.dist
 Description: Fit latent (hidden) Markov models on mixed categorical and continuous (timeseries)
    data, otherwise known as dependent mixture models

Modified: trunk/R/responseMVN.R
===================================================================
--- trunk/R/responseMVN.R	2010-02-24 10:58:56 UTC (rev 365)
+++ trunk/R/responseMVN.R	2010-02-24 14:49:50 UTC (rev 366)
@@ -1,4 +1,22 @@
-# Class 'MVNresponse' (multivariate normal response model)
+# Class 'MVNresponse' (multivariate normal response model)
+
+cov2par <- function(x) {
+	if(NROW(x) != NCOL(x)) stop("cov2par requires a square matrix") 
+	x[lower.tri(x,diag=TRUE)]
+}
+
+par2cov <- function(x) {
+	npar <- length(x)
+	dim <- (sqrt(8*npar + 1) - 1)/2
+	if(abs(dim - round(dim)) >= .Machine$double.eps^0.5) stop("number of parameters not suitable for par2cov")
+	cov <- matrix(0.0,ncol=dim,nrow=dim)
+	cov[lower.tri(cov,diag=TRUE)] <- x
+	cov[upper.tri(cov)] <- t(cov)[upper.tri(cov)]
+	cov
+}
+
+
+
 setClass("MVNresponse",
   representation(formula="formula"),
   contains="response"
@@ -10,7 +28,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 <- vech(cov.wt(x=fit$residuals,wt=w)$cov) else object at parameters$Sigma <- vech(cov(fit$residuals))
+		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
 	}
 )
@@ -69,13 +87,13 @@
 
 setMethod("logDens","MVNresponse",
 	function(object,...) {
-		dm_dmvnorm(x=object at y,mean=predict(object),sigma=xpnd(object at parameters$Sigma),log=TRUE,...)
+		dm_dmvnorm(x=object at y,mean=predict(object),sigma=par2cov(object at parameters$Sigma),log=TRUE,...)
 	}
 )
 
 setMethod("dens","MVNresponse",
 	function(object,log=FALSE,...) {
-		dm_dmvnorm(x=object at y,mean=predict(object),sigma=xpnd(object at parameters$Sigma),log=log,...)
+		dm_dmvnorm(x=object at y,mean=predict(object),sigma=par2cov(object at parameters$Sigma),log=log,...)
 	}
 )
 
@@ -96,11 +114,11 @@
 			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))
+		if(nrow(object at parameters$coefficients==1)) response <- mvrnorm(nt*nsim,mu=mu[1,],Sigma=par2cov(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))
+				response[i,] <- response <- mvrnorm(1,mu=mu[i,],Sigma=par2cov(object at parameters$Sigma))
 			}
 		}
 		return(response)
@@ -122,7 +140,7 @@
 		if(!is.matrix(y)) y <- matrix(y,ncol=1)
 		parameters <- list()
 		parameters$coefficients <- matrix(0.0,ncol=ncol(y),nrow=ncol(x))
-		parameters$Sigma <- vech(diag(ncol(y)))
+		parameters$Sigma <- cov2par(diag(ncol(y)))
 		npar <- length(unlist(parameters))		
 		if(is.null(fixed)) fixed <- as.logical(rep(0,npar))
 		if(!is.null(pstart)) {
@@ -142,7 +160,7 @@
 		cat("Coefficients: \n")
 		print(object at parameters$coefficients)
 		cat("Sigma: \n")
-		print(xpnd(object at parameters$Sigma))
+		print(par2cov(object at parameters$Sigma))
 	}
 )
 
@@ -181,4 +199,4 @@
 		)
 		return(pars)
 	}
-)
+)

Modified: trunk/man/depmix-internal.Rd
===================================================================
--- trunk/man/depmix-internal.Rd	2010-02-24 10:58:56 UTC (rev 365)
+++ trunk/man/depmix-internal.Rd	2010-02-24 14:49:50 UTC (rev 366)
@@ -18,6 +18,9 @@
 \alias{viterbi2}
 \alias{viterbi.fb}
 
+\alias{cov2par}
+\alias{par2cov}
+
 \alias{nlin}
 \alias{nlin,depmix.fitted-method}
 \alias{nlin,mix.fitted-method}