R/squeezeMVar.R

In betr: Identify differentially expressed genes in microarray time-course data

`squeezeMVar` <-
function (S, df, Lambda=NULL, nu=NULL) {
    n <- length(S)
    if (n == 0) 
        stop("Covariance matrix list is empty")
    if (n == 1) 
        return(list(varPrior = S[[1]], dfPrior = 0, varPost = S))
    k <- nrow(S[[1]])
    if (is.null(nu)) {
        Sdiag <- t(sapply(S, diag))    
        nu <- mean(sapply(1:k, function(x) squeezeVar(Sdiag[, x], df)$df.prior))
    }
    Savg <- matrix(apply(sapply(S, unlist), 1, mean), ncol=k)    
    nu.temp <- max (nu, k+6)
    if (is.null(Lambda)) { 
        if (is.infinite(nu)) {
            Lambda <- Savg 
        } else {
            Lambda <- (nu.temp-k-1) * Savg / nu.temp
        }
    }
    if (is.finite(nu)) {
        Stilde <- sapply(S, function(x) (df*x + nu*Lambda) / (df + nu), simplify=FALSE)
    } else {
        Stilde <- sapply(1:length(S), function(x) Lambda, simplify=FALSE)
    }
    Stilde <- lapply(Stilde, function(s) {
        rownames(s) <- colnames(s) <- NULL
        s})
    list(varPrior=Lambda, dfPrior=nu, varPost=Stilde)
}