[Depmix-commits] r185 - trunk/R

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

Author: maarten
Date: 2008-06-24 18:04:19 +0200 (Tue, 24 Jun 2008)
New Revision: 185

Modified:
   trunk/R/responseMVN.R
Log:
- added function to generate MVNresponse models
- fixed MVNnorm fit function (now also works without weights
- TODO: setpars function

Modified: trunk/R/responseMVN.R
===================================================================
--- trunk/R/responseMVN.R	2008-06-24 15:35:27 UTC (rev 184)
+++ trunk/R/responseMVN.R	2008-06-24 16:04:19 UTC (rev 185)
@@ -1,13 +1,17 @@
 # Class 'MVNresponse' (multivariate normal response model)
-setClass("MVNresponse",contains="response")
+setClass("MVNresponse",
+  representation(formula="formula"),
+  contains="response"
+)
 
 setMethod("fit","MVNresponse",
 	function(object,w) {
+    if(missing(w)) w <- NULL
 		pars <- object at parameters
-		fit <- lm.wfit(x=object at x,y=object at y,w=w)
+		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
-		object at parameters$Sigma <- cov.wt(x=fit$residuals,wt=w)["cov"]
-		object <- setpars(object,unlist(pars))
+		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))
 		object
 	}
 )
@@ -76,3 +80,46 @@
   }
 )
 
+MVNresponse <- function(formula,data,pstart=NULL,fixed=NULL,...) {
+	call <- match.call()
+	mf <- match.call(expand.dots = FALSE)
+	m <- match(c("formula", "data"), names(mf), 0)
+	mf <- mf[c(1, m)]
+	mf$drop.unused.levels <- TRUE
+	mf[[1]] <- as.name("model.frame")
+	mf <- eval(mf, parent.frame())
+	x <- model.matrix(attr(mf, "terms"),mf)
+	y <- model.response(mf)
+	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))
+	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))
+	   }
+  }
+	mod <- new("MVNresponse",formula=formula,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar)
+	mod
+}
+
+setMethod("show","MVNresponse",
+	function(object) {
+		cat("Multivariate Normal Model, formula: ",sep="")
+		print(object at formula)
+		cat("Coefficients: \n")
+		print(object at parameters$coefficients)
+		cat("Sigma: \n")
+    print(object at parameters$Sigma)
+	}
+)
+
+#TODO: setpars function