247 lines
8.7 KiB
R
247 lines
8.7 KiB
R
library(tensorPredictors)
|
|
suppressPackageStartupMessages(library(Rdimtools))
|
|
|
|
# Source utility function used in most simulations (extracted for convenience)
|
|
setwd("~/Work/tensorPredictors/sim/")
|
|
source("./sim_utils.R")
|
|
|
|
# Data set sample size in every simulation
|
|
sample.size <- 500L
|
|
# Nr. of per simulation replications
|
|
reps <- 10L
|
|
# number of observation/response axes (order of the tensors)
|
|
orders <- c(2L, 3L, 4L)
|
|
# auto correlation coefficient for the mode-wise auto scatter matrices `Omegas`
|
|
rhos <- seq(0, 0.8, by = 0.2)
|
|
|
|
|
|
setwd("~/Work/tensorPredictors/sim/")
|
|
base.name <- format(Sys.time(), "sim-tsir-%Y%m%dT%H%M")
|
|
# data sampling routine
|
|
sample.data <- function(sample.size, betas, Omegas) {
|
|
dimF <- mapply(ncol, betas)
|
|
|
|
# responce is a standard normal variable
|
|
y <- sort(rnorm(sample.size))
|
|
y.pow <- Reduce(function(a, b) outer(a, b, `+`), Map(seq, 0L, len = dimF))
|
|
F <- t(outer(y, as.vector(y.pow), `^`)) / as.vector(factorial(y.pow))
|
|
dim(F) <- c(dimF, sample.size)
|
|
|
|
matplot(mat(F, length(dim(F))), type = "l")
|
|
abline(h = 0, lty = "dashed")
|
|
|
|
matplot(y, scale(mat(F, length(dim(F))), scale = FALSE), type = "l")
|
|
abline(h = 0, lty = "dashed")
|
|
|
|
|
|
# sample predictors from tensor normal X | Y = y (last axis is sample axis)
|
|
sample.axis <- length(betas) + 1L
|
|
Sigmas <- Map(solve, Omegas)
|
|
mu_y <- mlm(F, Map(`%*%`, Sigmas, betas))
|
|
X <- mu_y + rtensornorm(sample.size, 0, Sigmas, sample.axis)
|
|
|
|
# Make `y` into a `Y` tensor with variable axis all of dim 1
|
|
Y <- array(y, dim = c(rep(1L, length(dimF)), sample.size))
|
|
|
|
list(X = X, F = F, Y = Y, sample.axis = sample.axis)
|
|
}
|
|
|
|
# Create a CSV logger to write simulation results to
|
|
log.file <- paste(base.name, "csv", sep = ".")
|
|
logger <- CSV.logger(
|
|
file.name = log.file,
|
|
header = c("rho", "order", "sample.size", "rep", "beta.version", outer(
|
|
"dist.subspace", c("gmlm", "gmlm.1d", "tsir", "sir"),
|
|
paste, sep = "."
|
|
))
|
|
)
|
|
|
|
for (order in orders) {
|
|
# setup problem dimensions given the tensor order
|
|
dimX <- rep(2L, order)
|
|
dimF <- rep(2L, order)
|
|
|
|
# as well as the projections onto rank 1 matrices
|
|
proj.betas <- Map(.projRank, rep(1L, order))
|
|
|
|
for (rho in rhos) {
|
|
# Scatter matrices (Omegas) given as autoregression structure `AR(rho)`
|
|
Omegas <- Map(function(p) rho^abs(outer(1:p, 1:p, `-`)), dimX)
|
|
|
|
# Version 1 betas (reductions)
|
|
beta.version <- 1L
|
|
betas <- Map(function(nr, nc) {
|
|
`[<-`(matrix(0, nr, nc), 1, , 1)
|
|
}, dimX, dimF)
|
|
B.true <- Reduce(kronecker, rev(betas))
|
|
# `(B.min %*% vec(X)) | Y = y` proportional to `(1 + y)^order`)
|
|
B.min.true <- B.true[, 1L, drop = FALSE]
|
|
|
|
# Version 1: repeated simulations
|
|
for (rep in seq_len(reps)) {
|
|
# Sample training data
|
|
c(X, F, Y, sample.axis) %<-% sample.data(sample.size, betas, Omegas)
|
|
|
|
# Fit models to provided data
|
|
fit.gmlm <- gmlm_tensor_normal(X, F, sample.axis = sample.axis, proj.betas = proj.betas)
|
|
fit.gmlm.y <- gmlm_tensor_normal(X, Y, sample.axis = sample.axis)
|
|
fit.tsir <- TSIR(X, drop(Y), d = rep(1L, order), sample.axis = sample.axis)
|
|
fit.sir <- do.sir(mat(X, sample.axis), drop(Y), ndim = 1L)
|
|
|
|
# Extract minimal reduction matrices from fitted models
|
|
B.gmlm <- qr.Q(qr(Reduce(kronecker, rev(fit.gmlm$betas))))[, 1L, drop = FALSE]
|
|
B.gmlm.y <- Reduce(kronecker, rev(fit.gmlm.y$betas))
|
|
B.tsir <- Reduce(kronecker, rev(fit.tsir))
|
|
B.sir <- fit.sir$projection
|
|
|
|
# Compute estimation to true minimal `B` distance
|
|
dist.subspace.gmlm <- dist.subspace(B.min.true, B.gmlm, normalize = TRUE)
|
|
dist.subspace.gmlm.y <- dist.subspace(B.min.true, B.gmlm.y, normalize = TRUE)
|
|
dist.subspace.tsir <- dist.subspace(B.min.true, B.tsir, normalize = TRUE)
|
|
dist.subspace.sir <- dist.subspace(B.min.true, B.sir, normalize = TRUE)
|
|
|
|
# Write to simulation log file (CSV file)
|
|
logger()
|
|
|
|
# and print progress
|
|
cat(sprintf("order %d, rho %.2f, version %d: rep: %d/%d\n",
|
|
order, rho, beta.version, rep, reps))
|
|
}
|
|
|
|
|
|
# Version 2 betas (reductions)
|
|
beta.version <- 2L
|
|
betas <- Map(function(nr, nc) {
|
|
tcrossprod((-1)^seq_len(nr), (-1)^seq_len(nc))
|
|
}, dimX, dimF)
|
|
B.true <- Reduce(kronecker, rev(betas))
|
|
# `(B.min %*% vec(X)) | Y = y` proportional to `(1 - y)^order`)
|
|
B.min.true <- B.true[, 1L, drop = FALSE]
|
|
|
|
# Version 2: repeated simulations (identical to Version 1)
|
|
for (rep in seq_len(reps)) {
|
|
# Sample training data
|
|
c(X, F, Y, sample.axis) %<-% sample.data(sample.size, betas, Omegas)
|
|
|
|
# Fit models to provided data
|
|
fit.gmlm <- gmlm_tensor_normal(X, F, sample.axis = sample.axis, proj.betas = proj.betas)
|
|
fit.gmlm.y <- gmlm_tensor_normal(X, Y, sample.axis = sample.axis)
|
|
fit.tsir <- TSIR(X, drop(Y), d = rep(1L, order), sample.axis = sample.axis)
|
|
fit.sir <- do.sir(mat(X, sample.axis), drop(Y), ndim = 1L)
|
|
|
|
# Extract minimal reduction matrices from fitted models
|
|
B.gmlm <- qr.Q(qr(Reduce(kronecker, rev(fit.gmlm$betas))))[, 1L, drop = FALSE]
|
|
B.gmlm.y <- Reduce(kronecker, rev(fit.gmlm.y$betas))
|
|
B.tsir <- Reduce(kronecker, rev(fit.tsir))
|
|
B.sir <- fit.sir$projection
|
|
|
|
# Compute estimation to true minimal `B` distance
|
|
dist.subspace.gmlm <- dist.subspace(B.min.true, B.gmlm, normalize = TRUE)
|
|
dist.subspace.gmlm.y <- dist.subspace(B.min.true, B.gmlm.y, normalize = TRUE)
|
|
dist.subspace.tsir <- dist.subspace(B.min.true, B.tsir, normalize = TRUE)
|
|
dist.subspace.sir <- dist.subspace(B.min.true, B.sir, normalize = TRUE)
|
|
|
|
# Write to simulation log file (CSV file)
|
|
logger()
|
|
|
|
# and print progress
|
|
cat(sprintf("order %d, rho %.2f, version %d: rep: %d/%d\n",
|
|
order, rho, beta.version, rep, reps))
|
|
}
|
|
}
|
|
}
|
|
|
|
### read simulation results and generate plots
|
|
if (!interactive()) { pdf(file = paste(base.name, "pdf", sep = ".")) }
|
|
|
|
# Read simulation results back from log file
|
|
sim <- read.csv(log.file)
|
|
|
|
# Source utility function used in most simulations (extracted for convenience)
|
|
source("./sim_utils.R")
|
|
|
|
# aggregate results
|
|
grouping <- sim[c("rho", "order", "beta.version")]
|
|
sim.dist <- sim[startsWith(names(sim), "dist")]
|
|
aggr <- aggregate(sim.dist, grouping, mean)
|
|
|
|
# reset (possible lost) config
|
|
orders <- sort(unique(aggr$order))
|
|
rhos <- sort(unique(aggr$rho))
|
|
|
|
# new grouping for the aggregates
|
|
layout(rbind(
|
|
matrix(seq_len(2 * length(orders)), ncol = 2),
|
|
2 * length(orders) + 1
|
|
), heights = c(rep(6L, length(orders)), 1L))
|
|
|
|
col.methods <- c(
|
|
gmlm = "#000000",
|
|
gmlm.y = "#FF0000",
|
|
tsir = "#009E73",
|
|
sir = "#999999"
|
|
)
|
|
|
|
for (group in split(aggr, aggr[c("order", "beta.version")])) {
|
|
order <- group$order[[1]]
|
|
beta.version <- group$beta.version[[1]]
|
|
|
|
rho <- group$rho
|
|
|
|
plot(range(rho), 0:1, main = sprintf("V%d, order %d", beta.version, order),
|
|
type = "n", bty = "n", axes = FALSE, xlab = expression(rho), ylab = "Subspace Distance")
|
|
axis(1, at = rho)
|
|
axis(2, at = seq(0, 1, by = 0.2))
|
|
with(group, {
|
|
lines(rho, dist.subspace.gmlm, col = col.methods["gmlm"], lwd = 3, type = "b", pch = 16)
|
|
lines(rho, dist.subspace.gmlm.y, col = col.methods["gmlm.y"], lwd = 3, type = "b", pch = 16)
|
|
lines(rho, dist.subspace.tsir, col = col.methods["tsir"], lwd = 2, type = "b", pch = 16)
|
|
lines(rho, dist.subspace.sir, col = col.methods["sir"], lwd = 2, type = "b", pch = 16)
|
|
})
|
|
if (order == 3L && beta.version == 2L) {
|
|
abline(v = 0.5, lty = "dotted", col = "black")
|
|
abline(h = group$dist.subspace.tsir[which(rho == 0.5)],
|
|
lty = "dotted", col = "black")
|
|
}
|
|
}
|
|
methods <- c("GMLM", "GMLM.y", "TSIR", "SIR")
|
|
restor.par <- par(
|
|
fig = c(0, 1, 0, 1),
|
|
oma = c(0, 0, 0, 0),
|
|
mar = c(1, 0, 0, 0),
|
|
new = TRUE
|
|
)
|
|
plot(0, type = "n", bty = "n", axes = FALSE, xlab = "", ylab = "")
|
|
legend("bottom", col = col.methods[tolower(methods)], legend = methods,
|
|
horiz = TRUE, lty = 1, bty = "n", lwd = c(3, 2, 2), pch = 16)
|
|
par(restor.par)
|
|
|
|
|
|
|
|
### Version 1
|
|
# 2'nd order
|
|
# 1 + 2 y + y^2
|
|
# (1 + y)^2
|
|
|
|
# 3'rd order
|
|
# 1 + 3 y + 3 y^2 + y^4
|
|
# (1 + y)^3
|
|
|
|
# 4'th order
|
|
# 1 + 4 y + 6 y^2 + 4 y^3 + y^4
|
|
# (1 + y)^4
|
|
|
|
|
|
### Version 2
|
|
# 2'nd order
|
|
# 1 - 2 y + y^2
|
|
# (1 - y)^2 = 1 - 2 y + y^2
|
|
|
|
# 3'rd order
|
|
# 1 - 3 y + 3 y^2 - y^3
|
|
# -(y - 1)^3
|
|
|
|
# 4'th order
|
|
# 1 - 4 y + 6 y^2 - 4 y^3 + y^4
|
|
# (y - 1)^4
|