1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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)
)
}
|
