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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
##############################################################
#' Fast simulation of r.v.s from a mixture of multivariate Student's t densities
#'
#' @param n number of random vectors to be simulated.
#' @param mu an (m x d) matrix, where m is the number of mixture components.
#' @param sigma as list of m covariance matrices (d x d) on for each mixture component.
#' Alternatively it can be a list of m cholesky decomposition of the covariance.
#' In that case \code{isChol} should be set to \code{TRUE}.
#' @param df a positive scalar representing the degrees of freedom. All the densities in the mixture have the same \code{df}.
#' @param w vector of length m, containing the weights of the mixture components.
#' @param ncores Number of cores used. The parallelization will take place only if OpenMP is supported.
#' @param isChol boolean set to true is \code{sigma} is the cholesky decomposition of the covariance matrix.
#' @param retInd when set to \code{TRUE} an attribute called "index" will be added to the output matrix of random variables.
#' This is a vector specifying to which mixture components each random vector belongs. \code{FALSE} by default.
#' @param A an (optional) numeric matrix of dimension (n x d), which will be used to store the output random variables.
#' It is useful when n and d are large and one wants to call \code{rmvn()} several times, without reallocating memory
#' for the whole matrix each time. NB: the element of \code{A} must be of class "numeric".
#' @param kpnames if \code{TRUE} the dimensions' names are preserved. That is, the i-th column of the output
#' has the same name as the i-th entry of \code{mu} or the i-th column of \code{sigma}.
#' \code{kpnames==FALSE} by default.
#' @return If \code{A==NULL} (default) the output is an (n x d) matrix where the i-th row is the i-th simulated vector.
#' If \code{A!=NULL} then the random vector are store in \code{A}, which is provided by the user, and the function
#' returns \code{NULL}. Notice that if \code{retInd==TRUE} an attribute called "index" will be added to A.
#' This is a vector specifying to which mixture components each random vector belongs.
#' @details There are many candidates for the multivariate generalization of Student's t-distribution, here we use
#' the parametrization described here \url{https://en.wikipedia.org/wiki/Multivariate_t-distribution}.
#'
#' Notice that this function does not use one of the Random Number Generators (RNGs) provided by R, but one
#' of the parallel cryptographic RNGs described in (Salmon et al., 2011). It is important to point out that this
#' RNG can safely be used in parallel, without risk of collisions between parallel sequence of random numbers.
#' The initialization of the RNG depends on R's seed, hence the \code{set.seed()} function can be used to
#' obtain reproducible results. Notice though that changing \code{ncores} causes most of the generated numbers
#' to be different even if R's seed is the same (see example below). NB: at the moment
#' the parallelization does not work properly on Solaris OS when \code{ncores>1}. Hence, \code{rmixt()} checks if the OS
#' is Solaris and, if this the case, it imposes \code{ncores==1}
#'
#' @author Matteo Fasiolo <matteo.fasiolo@@gmail.com>, C++ RNG engine by Thijs van den Berg <http://sitmo.com/>.
#' @references John K. Salmon, Mark A. Moraes, Ron O. Dror, and David E. Shaw (2011). Parallel Random Numbers: As Easy as 1, 2, 3.
#' D. E. Shaw Research, New York, NY 10036, USA.
#' @examples
#' # Create mixture of two components
#' df <- 6
#' mu <- matrix(rep(c(1, 2, 10, 20), 2), 2, 2, byrow = TRUE)
#' sigma <- list(diag(c(1, 10)), matrix(c(1, -0.9, -0.9, 1), 2, 2))
#' w <- c(0.1, 0.9)
#'
#' # Simulate
#' X <- rmixt(1e4, mu, sigma, df, w, retInd = TRUE)
#' plot(X, pch = '.', col = attr(X, "index"))
#'
#' # Simulate with fixed seed
#' set.seed(414)
#' rmixt(4, mu, sigma, df, w)
#'
#' set.seed(414)
#' rmixt(4, mu, sigma, df, w)
#'
#' set.seed(414)
#' rmixt(4, mu, sigma, df, w, ncores = 2) # r.v. generated on the second core are different
#'
#' ###### Here we create the matrix that will hold the simulated random variables upfront.
#' A <- matrix(NA, 4, 2)
#' class(A) <- "numeric" # This is important. We need the elements of A to be of class "numeric".
#'
#' set.seed(414)
#' rmixt(4, mu, sigma, df, w, ncores = 2, A = A) # This returns NULL ...
#' A # ... but the result is here
#'
#' @export rmixt
#'
rmixt <- function(n, mu, sigma, df, w, ncores = 1, isChol = FALSE, retInd = FALSE, A = NULL, kpnames = FALSE)
{
d <- ncol(mu)
m <- length(w)
if( !is.matrix(sigma) ) sigma <- as.matrix( sigma )
if( ncores > 1 && grepl('SunOS', Sys.info()['sysname']) ){
message("rmvn() cannot be used on multiple cores under Solaris. I am resetting \"ncores\" to 1.")
ncores <- 1
}
# Create output matrix
if( is.null(A) ) {
retMat <- TRUE # We return a matrix
A <- matrix(nrow = n, ncol = d)
class(A) <- "numeric"
} else {
retMat <- FALSE # We return NULL
if( class(A[1, 1]) != "numeric" ){
stop("class(A[1, 1]) != \"numeric\", to avoid this do class(A)<-\"numeric\".")
}
}
# Associate each sample with a mixture component
indV <- sample(0:(m-1), n, prob = w, replace = T)
.Call( "rmixtCpp",
n_ = n,
mu_ = mu,
sigma_ = sigma,
df_ = df,
indV_ = indV,
ncores_ = ncores,
isChol_ = isChol,
retInd_ = retInd,
A_ = A )
if( kpnames ) # (Optionally add dimensions names)
{
if( is.null(nam <- dimnames(mu)[[2L]]) ){ nam <- dimnames(sigma[[1L]])[[1L]] }
if( !is.null(nam) ){ dimnames(A) <- list(NULL, nam) }
}
# Return a matrix if no storage was provided and NULL if it was provided.
if( retMat ) {
return( A );
} else {
return( invisible(NULL) )
}
}
|
