# R/utilities.R In dmt: Dependency Modeling Toolkit

#### Documented in centerData

```# (C) 2008-2012 Leo Lahti and Olli-Pekka Huovilainen
# FreeBSD License (keep this notice)

# "Where a calculator on the ENIAC is equipped with 18,000 vacuum tubes
# and weighs 30 tons, computers in the future may have only 1,000
# vaccuum tubes and perhaps weigh 1.5 tons."
# - unknown, Popular Mechanics, March 1949

centerData <- function (X, rm.na = TRUE, meanvalue = NULL) {

# Shift data matrix (columns) to zero, or given 'meanvalue'

if (!rm.na) {
xcenter <- colMeans(X)
X2 <- X - rep(xcenter, rep.int(nrow(X), ncol(X)))
} else {
X2 <- array(NA, dim = c(nrow(X), ncol(X)), dimnames = dimnames(X))
for (i in 1:ncol(X)) {
x <- X[,i]
nainds <- is.na(x)
xmean <- mean(x[!nainds])
X2[!nainds,i] <- x[!nainds] - xmean
}
dimnames(X2) <- dimnames(X)
}

if (!is.null(meanvalue)) {
# Shift the data so that mean gets a specified value
X2 <- X2 + meanvalue
}

X2
}

initialize2 <- function (X, Y, zDim = NULL, marginalCovariances) {

zDim <- ifelse(is.null(zDim), min(nrow(X), nrow(Y)), zDim)

Nsamples <- ncol(X)
Dim <- list(X = nrow(X), Y = nrow(Y), Z = zDim)
nullmat  <- matrix(0, nrow = Dim\$X, ncol = Dim\$Y)

if (marginalCovariances == "isotropic") {
# Scalar values
phi <- list()
phi\$X <- diag(var(as.vector(X)), Dim\$X)
phi\$Y <- diag(var(as.vector(Y)), Dim\$Y)
phi\$total <- diag(c(diag(phi\$X), diag(phi\$Y)))
#nullmat <- 0
} else if (marginalCovariances == "identical isotropic") {
# Scalar values phix = phiy
phi <- list()
phi.est <- var(c(as.vector(X), as.vector(Y)))
phi\$X <- diag(phi.est, Dim\$X)
phi\$Y <- diag(phi.est, Dim\$Y)
phi\$total <- diag(c(diag(phi\$X), diag(phi\$Y)))
#nullmat <- 0
} else {
# diagonal matrices
# initialize with scalar diagonal noise on the marginals (shared by all features)
phi <- list(X = diag(var(as.vector(X)), Dim\$X),
Y = diag(var(as.vector(Y)), Dim\$Y))
phi\$total <- rbind(cbind(phi\$X,nullmat), cbind(t(nullmat), phi\$Y))

}

# FIXME: if phi\$Y is scalar (as in segmented/mir case) we can speed up here. Do later.

phi.inv  <- list()
phi.inv\$X <- solve(phi\$X)
phi.inv\$Y <- solve(phi\$Y)
#phi.inv\$total <- rbind(cbind(phi.inv\$X, nullmat), cbind(t(nullmat), phi.inv\$Y))

Dcov <- list()
Dcov\$X <- cov(t(X), use = "pairwise.complete.obs")
Dcov\$Y <- cov(t(Y), use = "pairwise.complete.obs")
Dcov\$xy <- cov(t(X), t(Y), use = "pairwise.complete.obs")
Dcov\$yx <- t(Dcov\$xy)
Dcov\$total <- rbind(cbind(Dcov\$X, Dcov\$xy), cbind(Dcov\$yx, Dcov\$Y))
if (nrow(X) == nrow(Y)) {
Dcov\$sum <- Dcov\$X + Dcov\$Y + Dcov\$xy + Dcov\$yx
#Dcov\$sum <- cov(t(X + Y), use = "pairwise.complete.obs")
}

# It is possible that covariances calculated with pairwise complete
# observations are not positive semi-definite.
# Check this. If not pos.sem.def, then replace with the closest
# pos. semidefinite matrix
if (any(eigen(Dcov\$total)\$values < 0)) {
#message("Covariance approximation used to avoid numerical instability.")
Dcov\$X  <- as.matrix(nearPD(Dcov\$X)\$mat)
Dcov\$Y  <- as.matrix(nearPD(Dcov\$Y)\$mat)
Dcov\$total <- rbind(cbind(Dcov\$X, Dcov\$xy), cbind(Dcov\$yx, Dcov\$Y))
if (nrow(X) == nrow(Y)) {
Dcov\$sum <- Dcov\$X + Dcov\$Y + Dcov\$xy + Dcov\$yx
}
}

# Initialize W's
W <- list()
W\$X <- as.matrix(eigen(Dcov\$X)\$vectors[, 1:Dim\$Z])
W\$Y <- as.matrix(eigen(Dcov\$Y)\$vectors[, 1:Dim\$Z])
W\$total <- rbind(W\$X, W\$Y)

list(phi = phi, W = W, Dcov = Dcov, Dim = Dim, nullmat = nullmat, Nsamples = Nsamples)

}
```

## Try the dmt package in your browser

Any scripts or data that you put into this service are public.

dmt documentation built on May 29, 2017, 11 a.m.