add: special case for vector valued F in gmlm_tensor_normal()
fix: typos
This commit is contained in:
parent
6cbb5029d4
commit
57eb279e6a
|
|
@ -12,4 +12,4 @@ License: GPL (>= 2)
|
|||
Imports: Rcpp (>= 1.0.8)
|
||||
LinkingTo: Rcpp
|
||||
SystemRequirements: c++17
|
||||
RoxygenNote: 7.2.0
|
||||
RoxygenNote: 7.3.2
|
||||
|
|
|
|||
285
sim/sim-tsir.R
285
sim/sim-tsir.R
|
|
@ -1,246 +1,127 @@
|
|||
library(tensorPredictors)
|
||||
suppressPackageStartupMessages(library(Rdimtools))
|
||||
|
||||
setwd("~/Work/tensorPredictors/sim/")
|
||||
base.name <- format(Sys.time(), "sim_tsir-%Y%m%dT%H%M")
|
||||
|
||||
# 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)
|
||||
# Set PRNG seed for reproducability
|
||||
# Sequence 'T', 'S', 'I', 'R' in ASCII is 84, 83, 73, 82.
|
||||
set.seed(84837382L, "Mersenne-Twister", "Inversion", "Rejection")
|
||||
|
||||
### Simulation configuration
|
||||
reps <- 100 # number of simulation replications
|
||||
sample.sizes <- c(100, 200, 300, 500, 750) # sample sizes `n`
|
||||
signal.noise.ratios <- 2^(-3:4) # Signal to Noise Ratios (from 50/50 to very high)
|
||||
dimX <- c(2, 3, 5) # predictor `X` dimension
|
||||
dimF <- rep(2, length(dimX)) # "function" `F(y)` of responce `y` dimension
|
||||
|
||||
# setup true model parameters
|
||||
eta1 <- 0
|
||||
# rank 1 betas
|
||||
betas <- Map(function(nr, nc) {
|
||||
tcrossprod((-1)^seq_len(nr), (-1)^seq_len(nc))
|
||||
}, dimX, dimF)
|
||||
# True (minimal) reduction matrix
|
||||
B.true <- Reduce(kronecker, rev(betas))[, 1L, drop = FALSE]
|
||||
# GMLM second moment parameters (mode-wise precition matrices)
|
||||
Omegas <- Map(function(pj) 0.5^abs(outer(1:pj, 1:pj, `-`)), dimX) # AR(0.5)
|
||||
# True (minimal) Gamma (Projection Direction)
|
||||
Gamma.true <- Reduce(kronecker, rev(Map(solve, Omegas, betas)))[, 1L, drop = FALSE]
|
||||
# true (full) covariance matrix
|
||||
covX.true <- Reduce(kronecker, rev(Map(solve, Omegas)))
|
||||
|
||||
# define projections onto rank 1 matrices for betas
|
||||
proj.betas <- Map(.projRank, rep(1L, length(betas)))
|
||||
|
||||
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)
|
||||
|
||||
sample.data <- function(sample.size, eta1, betas, Omegas, snr) {
|
||||
# responce is a standard normal variable
|
||||
y <- sort(rnorm(sample.size))
|
||||
y <- rnorm(sample.size)
|
||||
# F(y) is a tensor of monomials
|
||||
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))
|
||||
F <- t(outer(y, as.vector(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)
|
||||
Deltas <- Map(solve, Omegas) # normal covariances
|
||||
mu_y <- mlm(mlm(F, betas) + as.vector(eta1), Deltas) # conditional mean
|
||||
noise <- rtensornorm(sample.size, 0, Deltas, sample.axis) # error term
|
||||
# scale noise to given signal to noise ratio
|
||||
snr.est <- sd(mu_y) / sd(noise)
|
||||
noise <- (snr.est / snr) * noise
|
||||
X <- mu_y + noise # responses X
|
||||
|
||||
# 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)
|
||||
list(X = X, F = F, y = y, sample.axis = sample.axis)
|
||||
}
|
||||
|
||||
# Create a CSV logger to write simulation results to
|
||||
# Create a CSV logger to save simulation results
|
||||
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 = "."
|
||||
))
|
||||
header = c(
|
||||
"snr", "sample.size", "rep",
|
||||
"dist.subspace.gmlm", "dist.subspace.tsir", "dist.subspace.partial", "dist.subspace.gamma"
|
||||
)
|
||||
)
|
||||
|
||||
for (order in orders) {
|
||||
# setup problem dimensions given the tensor order
|
||||
dimX <- rep(2L, order)
|
||||
dimF <- rep(2L, order)
|
||||
# different Signal to Noise Ratios
|
||||
for (snr in signal.noise.ratios) {
|
||||
|
||||
# as well as the projections onto rank 1 matrices
|
||||
proj.betas <- Map(.projRank, rep(1L, order))
|
||||
# simulation for multiple data set sizes
|
||||
for (sample.size in sample.sizes) {
|
||||
|
||||
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
|
||||
# simulation replications
|
||||
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)
|
||||
# sample a data set
|
||||
c(X, F, y, sample.axis) %<-% sample.data(sample.size, eta1, betas, Omegas, snr)
|
||||
|
||||
# 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
|
||||
# call GMLM and TSIR
|
||||
fit.gmlm <- gmlm_tensor_normal(X, F, sample.axis = sample.axis, proj.betas = proj.betas)
|
||||
fit.tsir <- TSIR(X, y, c(1L, 1L, 1L), sample.axis = sample.axis)
|
||||
|
||||
# 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)
|
||||
# GMLM, TSIR reduction estimates and TSIR (internal) projections
|
||||
B.gmlm <- Reduce(kronecker, Map(function(beta) qr.Q(qr(beta))[, 1L, drop = FALSE], rev(fit.gmlm$betas)))
|
||||
B.tsir <- Reduce(kronecker, rev(fit.tsir))
|
||||
Gamma <- Reduce(kronecker, rev(attr(fit.tsir, "Gammas")))
|
||||
|
||||
# Subspace distances to true minimal reduction
|
||||
dist.subspace.gmlm <- dist.subspace(B.true, B.gmlm, normalize = TRUE)
|
||||
dist.subspace.tsir <- dist.subspace(B.true, B.tsir, normalize = TRUE)
|
||||
dist.subspace.partial <- dist.subspace(B.true, solve(covX.true, Gamma), normalize = TRUE)
|
||||
dist.subspace.gamma <- dist.subspace(Gamma.true, Gamma, 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))
|
||||
cat(sprintf("SNR %.2f, sample size %d: rep: %d/%d\n", snr, sample.size, rep, reps))
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
### read simulation results and generate plots
|
||||
|
||||
### read simulation results generate plots
|
||||
if (!interactive()) { pdf(file = paste(base.name, "pdf", sep = ".")) }
|
||||
|
||||
# Read simulation results back from log file
|
||||
# Read siulation results from log file
|
||||
sim <- read.csv(log.file)
|
||||
|
||||
# Source utility function used in most simulations (extracted for convenience)
|
||||
source("./sim_utils.R")
|
||||
# reset the correlation configuration parameter
|
||||
signal.noise.ratios <- sort(unique(sim$snr))
|
||||
|
||||
# aggregate results
|
||||
grouping <- sim[c("rho", "order", "beta.version")]
|
||||
sim.dist <- sim[startsWith(names(sim), "dist")]
|
||||
aggr <- aggregate(sim.dist, grouping, mean)
|
||||
# build plot layout for every `snr` param
|
||||
ncols <- ceiling(sqrt(length(signal.noise.ratios)))
|
||||
nrows <- ceiling(length(signal.noise.ratios) / ncols)
|
||||
par(mfrow = c(nrows, ncols))
|
||||
|
||||
# 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")
|
||||
}
|
||||
# One plot for every Singal to Noise Ratio
|
||||
for (.snr in signal.noise.ratios) {
|
||||
plot.sim(subset(sim, snr == .snr), "dist.subspace",
|
||||
main = sprintf("Signal to Noise Ratio: %.3f", .snr), xlab = "Sample Size", ylab = "Subspace Distance",
|
||||
cols = c(gmlm = "black", tsir = "#009E73", partial = "orange", gamma = "skyblue")
|
||||
)
|
||||
}
|
||||
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
|
||||
|
|
|
|||
|
|
@ -1,4 +1,6 @@
|
|||
library(tensorPredictors)
|
||||
library(rTensor)
|
||||
|
||||
# library(RGCCA)
|
||||
### Load modified version which _does not_ create a clusster in case of
|
||||
### `n_cores == 1` allowing huge speed improvements! (at least on Ubuntu 22.04 LTS)
|
||||
|
|
|
|||
|
|
@ -25,17 +25,15 @@ dimX <- c(2, 3, 5) # predictor `X` dimension
|
|||
dimF <- rep(2, length(dimX)) # "function" `F(y)` of responce `y` dimension
|
||||
|
||||
# setup true model parameters (rank 1 betas)
|
||||
betas <- list(
|
||||
`[<-`(matrix(0, dimX[1], dimF[1]), 1, , c(1, 1)),
|
||||
`[<-`(matrix(0, dimX[2], dimF[2]), 2, , c(1, 0)),
|
||||
`[<-`(matrix(0, dimX[3], dimF[3]), 3, , c(1, 2))
|
||||
)
|
||||
betas <- Map(function(nr, nc) {
|
||||
tcrossprod((-1)^seq_len(nr), (-1)^seq_len(nc))
|
||||
}, dimX, dimF)
|
||||
# invisible(Map(print.table, betas, zero.print = "."))
|
||||
Omegas <- Map(function(pj) 0.5^abs(outer(1:pj, 1:pj, `-`)), dimX) # AR(0.5)
|
||||
eta1 <- 0
|
||||
|
||||
# True (full) reduction matrix to compair against
|
||||
B.true <- as.matrix(as.numeric(15 == (1:30)))
|
||||
B.true <- Reduce(kronecker, rev(betas))[, 1L, drop = FALSE]
|
||||
|
||||
# define projections onto rank 1 matrices for betas
|
||||
proj.betas <- Map(.projRank, rep(1L, length(betas)))
|
||||
|
|
@ -70,46 +68,6 @@ logger <- CSV.logger(
|
|||
))
|
||||
)
|
||||
|
||||
# # true reduction (1D, select first component)
|
||||
# B.true <- as.matrix(as.numeric(1 == (1:30)))
|
||||
|
||||
|
||||
|
||||
# B.true <- Reduce(`%x%`, rev(betas)) #[, 1L, drop = FALSE]
|
||||
# print.table(B.true, zero.print = ".")
|
||||
|
||||
# print.table(Reduce(kronecker, rev(betas)), zero.print = ".")
|
||||
# mu1 <- function(y) 1 + 3 * y + y^2
|
||||
|
||||
# matX <- mat(X, 1:3)
|
||||
# x1 <- matX[1, ]
|
||||
# x2 <- matX[2, ]
|
||||
|
||||
# plot(mu1(y), x1, col = cut(y, 10L))
|
||||
# plot(y, x1, col = cut(y, 10L))
|
||||
|
||||
# plot(mu1(y), x2, col = cut(y, 10L))
|
||||
# plot(y, matX[1, ], col = cut(y, 10L))
|
||||
# plot(y, matX[10, ], col = cut(y, 10L))
|
||||
|
||||
# c(X, F, y, sample.axis) %<-% sample.data(sample.size, eta1, betas, Omegas)
|
||||
|
||||
|
||||
# colnames(B.true) <- c("1", "y", "y", "y^2", "y", "y^2", "y^2", "y^3")
|
||||
|
||||
# mu_y <- function(y) {
|
||||
# B.min <- cbind(
|
||||
# B.true[, 1],
|
||||
# rowSums(B.true[, c(2:3, 5)]),
|
||||
# rowSums(B.true[, c(4, 6:7)]),
|
||||
# B.true[, 8]
|
||||
# )
|
||||
|
||||
# B.min %*% rbind(1, y, y^2, y^3)
|
||||
# }
|
||||
|
||||
# plot((mat(X, 4) %*% B.min)[, 2], y)
|
||||
|
||||
### for each sample size
|
||||
for (sample.size in sample.sizes) {
|
||||
# repeate every simulation
|
||||
|
|
@ -139,7 +97,7 @@ for (sample.size in sample.sizes) {
|
|||
F1 <- cbind(y, y^2, y^3)
|
||||
time.mgcca <- system.time(
|
||||
fit.mgcca <- mgcca(
|
||||
list(X1, F1), # `drop` removes 1D axis
|
||||
list(X1, F1),
|
||||
quiet = TRUE,
|
||||
scheme = "factorial",
|
||||
ncomp = c(1L, 1L)
|
||||
|
|
|
|||
|
|
@ -37,7 +37,7 @@ betas <- Map(diag, 1, dimX, dimF)
|
|||
# Omegas <- Map(function(dim) `diag<-`(matrix(1, dim, dim), 0), dimX)
|
||||
Omegas <- list(
|
||||
`diag<-`(matrix(0.5, dimX[1], dimX[1]), 0),
|
||||
diag(dimX[2])
|
||||
matrix(c(0, 1, 0, 1, 0, 1, 0, 1, 0), 3, 3)
|
||||
)
|
||||
|
||||
# compute true (full) model parameters to compair estimates against
|
||||
|
|
|
|||
|
|
@ -142,13 +142,11 @@ names(col.methods) <- methods
|
|||
|
||||
|
||||
# Comparison plot of one measure for a simulation
|
||||
plot.sim <- function(sim, measure.name, ..., ylim = c(0, 1)) {
|
||||
plot.sim <- function(sim, measure.name, ..., ylim = c(0, 1), cols = col.methods) {
|
||||
par.default <- par(pch = 16, bty = "n", lty = "solid", lwd = 1.5)
|
||||
|
||||
# # Set colors for every method
|
||||
# methods <- c("gmlm", "pca", "hopca", "tsir", "mgcca", "lpca", "clpca", "tnormal")
|
||||
# col.methods <- palette.colors(n = length(methods), palette = "Okabe-Ito", recycle = FALSE)
|
||||
# names(col.methods) <- methods
|
||||
# extract method names from color specification
|
||||
methods <- names(cols)
|
||||
|
||||
# Remain sample size grouping variable to avoid conflicts
|
||||
aggr.mean <- aggregate(sim, list(sampleSize = sim$sample.size), mean)
|
||||
|
|
@ -166,7 +164,7 @@ plot.sim <- function(sim, measure.name, ..., ylim = c(0, 1)) {
|
|||
min <- aggr.min[aggr.sd$sampleSize == sampleSize, dist.name]
|
||||
max <- aggr.max[aggr.sd$sampleSize == sampleSize, dist.name]
|
||||
method <- tail(strsplit(dist.name, ".", fixed = TRUE)[[1]], 1)
|
||||
col <- col.methods[method]
|
||||
col <- cols[method]
|
||||
lines(sampleSize, mean, type = "o", col = col, lty = 1, lwd = 2 + (method == "gmlm"))
|
||||
lines(sampleSize, mean + sd, col = col, lty = 2, lwd = 0.8)
|
||||
lines(sampleSize, mean - sd, col = col, lty = 2, lwd = 0.8)
|
||||
|
|
@ -175,7 +173,7 @@ plot.sim <- function(sim, measure.name, ..., ylim = c(0, 1)) {
|
|||
lines(sampleSize, max, col = col, lty = 3, lwd = 0.6)
|
||||
}
|
||||
|
||||
legend("topright", col = col.methods, lty = 1, legend = names(col.methods),
|
||||
legend("topright", col = cols, lty = 1, legend = names(cols),
|
||||
bty = "n", lwd = par("lwd"), pch = par("pch"))
|
||||
})
|
||||
|
||||
|
|
|
|||
|
|
@ -7,6 +7,11 @@ gmlm_tensor_normal <- function(X, F, sample.axis = length(dim(X)),
|
|||
max.iter = 100L, proj.betas = NULL, proj.Omegas = NULL, logger = NULL,
|
||||
cond.threshold = 25, eps = 1e-6
|
||||
) {
|
||||
# Special case for `F` being vector valued, add dims of size 1
|
||||
if (is.null(dim(F))) {
|
||||
dim(F) <- ifelse(seq_along(dim(X)) == sample.axis, length(F), 1L)
|
||||
}
|
||||
|
||||
# rearrange `X`, `F` such that the last axis enumerates observations
|
||||
if (!missing(sample.axis)) {
|
||||
axis.perm <- c(seq_along(dim(X))[-sample.axis], sample.axis)
|
||||
|
|
|
|||
|
|
@ -35,7 +35,7 @@ mcov <- function(X, sample.axis = length(dim(X)), center = TRUE) {
|
|||
X <- X - c(rowMeans(X, dims = r))
|
||||
}
|
||||
|
||||
# estimes (unscaled) covariances for each mode
|
||||
# unscaled covariances for each mode
|
||||
Sigmas <- .mapply(mcrossprod, list(mode = seq_len(r)), MoreArgs = list(X))
|
||||
# scale by per mode "sample" size
|
||||
Sigmas <- .mapply(`*`, list(Sigmas, p / prod(dim(X))), NULL)
|
||||
|
|
|
|||
|
|
@ -20,6 +20,9 @@ extern SEXP R_mcrossprod(SEXP A, SEXP B, SEXP mode);
|
|||
/* Symmetric Tensor Mode Crossproduct `A_(m) A_(m)^T` */
|
||||
extern SEXP R_mcrossprod_sym(SEXP A, SEXP mode);
|
||||
|
||||
/* Merge matrix-matrix multiplication over 3rd axis of 3D arrays */
|
||||
extern SEXP R_merge_matmul(SEXP A, SEXP B);
|
||||
|
||||
// /* Higher Order PCA */
|
||||
// extern SEXP hopca(SEXP X);
|
||||
|
||||
|
|
@ -73,6 +76,7 @@ static const R_CallMethodDef CallEntries[] = {
|
|||
{"C_mtvk", (DL_FUNC) &mtvk, 2},
|
||||
{"C_mcrossprod", (DL_FUNC) &R_mcrossprod, 3},
|
||||
{"C_mcrossprod_sym", (DL_FUNC) &R_mcrossprod_sym, 2},
|
||||
{"C_merge_matmul", (DL_FUNC) &R_merge_matmul, 2},
|
||||
{"C_svd", (DL_FUNC) &R_svd, 1},
|
||||
{"C_solve", (DL_FUNC) &R_solve, 2},
|
||||
{"C_det", (DL_FUNC) &R_det, 1},
|
||||
|
|
|
|||
|
|
@ -419,7 +419,7 @@ extern SEXP R_ising_m2(
|
|||
} else {
|
||||
// Exact method (ignore other arguments), only validate dimension
|
||||
if (25 < dim) {
|
||||
Rf_error("Dimension '%d' too big for exact method (max 24)", dim);
|
||||
Rf_error("Dimension '%ld' too big for exact method (max 24)", dim);
|
||||
}
|
||||
|
||||
// and call the exact method
|
||||
|
|
|
|||
Loading…
Reference in New Issue