[Pomp-commits] r148 - in pkg: . R

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

Author: kingaa
Date: 2009-06-17 18:08:24 +0200 (Wed, 17 Jun 2009)
New Revision: 148

Removed:
   pkg/R/mvnorm.R
Modified:
   pkg/DESCRIPTION
   pkg/NAMESPACE
   pkg/R/nlf-guts.R
Log:
make use of mvtnorm::dmvnorm now instead of our own version

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2009-06-17 15:46:35 UTC (rev 147)
+++ pkg/DESCRIPTION	2009-06-17 16:08:24 UTC (rev 148)
@@ -1,11 +1,11 @@
 Package: pomp
 Type: Package
 Title: Statistical inference for partially observed Markov processes
-Version: 0.24-5
-Date: 2009-06-03
+Version: 0.24-6
+Date: 2009-06-17
 Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
 Maintainer: Aaron A. King <kingaa at umich.edu>
 Description: Inference methods for partially-observed Markov processes
-Depends: R(>= 2.7.1), stats, methods, graphics, deSolve, subplex
+Depends: R(>= 2.8.1), stats, methods, graphics, deSolve, subplex, mvtnorm
 License: GPL (>= 2)
 LazyLoad: true

Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE	2009-06-17 15:46:35 UTC (rev 147)
+++ pkg/NAMESPACE	2009-06-17 16:08:24 UTC (rev 148)
@@ -17,6 +17,7 @@
 
 importFrom(graphics,plot)		
 importFrom(stats,simulate,time,coef,logLik)
+importFrom(mvtnorm,dmvnorm)
 
 exportClasses('pomp','mif')
 

Deleted: pkg/R/mvnorm.R
===================================================================
--- pkg/R/mvnorm.R	2009-06-17 15:46:35 UTC (rev 147)
+++ pkg/R/mvnorm.R	2009-06-17 16:08:24 UTC (rev 148)
@@ -1,72 +0,0 @@
-## multivariate normal density
-## stolen from package 'mvtnorm' version 0.9-7
-
-
-rmvnorm <- function (n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)),
-                     method=c("eigen", "svd", "chol"))
-{    
-    if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps))) {
-        stop("sigma must be a symmetric matrix")
-    }
-    if (length(mean) != nrow(sigma)) {
-        stop("mean and sigma have non-conforming size")
-    }
-    sigma1 <- sigma
-    dimnames(sigma1) <- NULL
-    if(!isTRUE(all.equal(sigma1, t(sigma1)))){
-        warning("sigma is numerically not symmetric")
-    }
-
-    method <- match.arg(method)
-    
-    if(method == "eigen"){
-        ev <- eigen(sigma, symmetric = TRUE)
-        if (!all(ev$values >= -sqrt(.Machine$double.eps) * abs(ev$values[1]))){
-            warning("sigma is numerically not positive definite")
-        }    
-        retval <- ev$vectors %*%  diag(sqrt(ev$values),length(ev$values)) %*% t(ev$vectors)
-    }
-    else if(method == "svd"){
-        sigsvd <- svd(sigma)
-        if (!all(sigsvd$d >= -sqrt(.Machine$double.eps) * abs(sigsvd$d[1]))){
-            warning("sigma is numerically not positive definite")
-        }    
-        retval <- t(sigsvd$v %*% (t(sigsvd$u) * sqrt(sigsvd$d)))
-    }    
-    else if(method == "chol"){
-        retval <- chol(sigma, pivot = TRUE)
-        o <- order(attr(retval, "pivot"))
-        retval <- retval[,o]
-    }
-    
-    retval <- matrix(rnorm(n * ncol(sigma)), nrow = n) %*%  retval
-    retval <- sweep(retval, 2, mean, "+")
-    retval
-}
-
-dmvnorm <- function (x, mean, sigma, log=FALSE)
-{
-    if (is.vector(x)) {
-        x <- matrix(x, ncol = length(x))
-    }
-    if (missing(mean)) {
-        mean <- rep(0, length = ncol(x))
-    }
-    if (missing(sigma)) {
-        sigma <- diag(ncol(x))
-    }
-    if (NCOL(x) != NCOL(sigma)) {
-        stop("x and sigma have non-conforming size")
-    }
-    if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps))) {
-        stop("sigma must be a symmetric matrix")
-    }
-    if (length(mean) != NROW(sigma)) {
-        stop("mean and sigma have non-conforming size")
-    }
-    distval <- mahalanobis(x, center = mean, cov = sigma)
-    logdet <- sum(log(eigen(sigma, symmetric=TRUE, only.values=TRUE)$values))
-    logretval <- -(ncol(x)*log(2*pi) + logdet + distval)/2
-    if(log) return(logretval)
-    exp(logretval)
-}

Modified: pkg/R/nlf-guts.R
===================================================================
--- pkg/R/nlf-guts.R	2009-06-17 15:46:35 UTC (rev 147)
+++ pkg/R/nlf-guts.R	2009-06-17 16:08:24 UTC (rev 148)
@@ -185,7 +185,7 @@
     LQL <- dnorm(prediction.errors[,1],mean=0,sd=sigma.model,log=TRUE)
   } else {
     sigma.model <- cov(model.residuals)
-    LQL <- dmvnorm(prediction.errors,sigma=sigma.model,log=TRUE) ## NOTE: This could be improved using GLS.
+    LQL <- mvtnorm::dmvnorm(prediction.errors,sigma=sigma.model,log=TRUE) ## NOTE: This could be improved using GLS.
   }  
   
   LQL