Start R package for the RSC simulation and analysis

1 files changed, 87 insertions, 0 deletions

diff --git a/R/simulate_data.R b/R/simulate_data.R
new file mode 100644
index 0000000..2338805
--- /dev/null
+++ b/R/simulate_data.R
@@ -0,0 +1,87 @@
+#' Simulate Multivariate Normal Data
+#'
+#' Given a set of parameters controlling the properties
+#' of the true covariance matrix, simulate data
+#' with the intention of evaluating the performance
+#' of covariance or correlation matrix estimators.
+#'
+#' @param n Number of observations to generate.
+#' @param p Number of features to simulate.
+#' @param rho The AR(1) decay parameter (default 0.8). Setting rho = 0 uses the identity matrix as the true covariance (0 <= rho < 1).
+#' @param contamination The proportions of rows to contaminate as outlier observations (default 0, 0 <= contamination <= 1).
+#' @param block_fraction The fraction of the covariance matrix to form as AR(1) structure (default 0.5, 0 < block_fraction <= 1).
+#' @param t_df The degrees of freedom to use for scaling observations into heavy tails (default Inf, 0 <= t_df < 2).
+#'
+#' @return A list with X, S and a vector identifying the contaminated rows.
+#'
+#' @examples
+#' X_obj <- generate_data(n = 10, p = 10)
+#' X <- X_obj$X
+#' S_true <- X_obj$S
+#'
+#' @importFrom MASS mvrnorm
+#'
+#' @export
+simulate_data <- function(
+	n,
+	p,
+	rho = 0.8,
+	contamination = 0,
+	block_fraction = 0.5,
+	t_df = Inf,
+	exact_contamination = FALSE
+){
+	if(any(
+		rho < 0,
+		rho >= 1,
+		contamination < 0,
+		contamination > 1,
+		block_fraction <= 0,
+		block_fraction > 1,
+		t_df <= 2
+	)){
+		stop('Incompatible settings')
+	}
+
+	# True correlation matrix
+	Sigma <- diag(p)
+	if(rho > 0){
+		b <- floor(block_fraction * p)
+		if(b >= 2){
+			idx <- seq_len(b)
+			Sigma[idx, idx] <- rho ^ abs(outer(idx, idx, `-`))
+		}
+	}
+
+	# Regular clean observations
+	Z <- MASS::mvrnorm(n = n, mu = numeric(p), Sigma = Sigma)
+	if(is.finite(t_df)){
+		w <- rchisq(n, df = t_df)
+		X <- Z / sqrt(w / t_df)
+	}else{
+		X <- Z
+	}
+
+	# Contamination - by row
+	outlier_rows <- integer(0)
+	if(contamination > 0){
+		if(exact_contamination){
+			outlier_rows <- sample.int(n, ceiling(contamination * n))
+		}else{
+			outlier_rows <- which(runif(n) < contamination)
+		}
+
+		if(length(outlier_rows) > 0){
+			X[outlier_rows, ] <- matrix(
+				rchisq(length(outlier_rows) * p, df = 1),
+				nrow = length(outlier_rows),
+				ncol = p
+			)
+		}
+	}
+	list(
+		X = X,
+		S = Sigma,
+		contaminated = sort(outlier_rows)
+	)
+}