Fast density computation for mixture of multivariate normal distributions.

Fast density computation for mixture of multivariate normal distributions.

dmixn(X, mu, sigma, w, log = FALSE, ncores = 1, isChol = FALSE,
  A = NULL)

Arguments

X	matrix n by d where each row is a d dimensional random vector. Alternatively X can be a d-dimensional vector.
mu	an (m x d) matrix, where m is the number of mixture components.
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 isChol should be set to TRUE.
w	vector of length m, containing the weights of the mixture components.
log	boolean set to true the logarithm of the pdf is required.
ncores	Number of cores used. The parallelization will take place only if OpenMP is supported.
isChol	boolean set to true is sigma is the cholesky decomposition of the covariance matrix.
A	an (optional) numeric matrix of dimension (m x d), which will be used to store the evaluations of each mixture density over each mixture component. It is useful when m and n are large and one wants to call dmixt() several times, without reallocating memory for the whole matrix each time. NB1: A will be modified, not copied! NB2: the element of A must be of class "numeric".

Value

A vector of length n where the i-the entry contains the pdf of the i-th random vector (i.e. the i-th row of X).

Details

NB: at the moment the parallelization does not work properly on Solaris OS when ncores>1. Hence, dmixt() checks if the OS is Solaris and, if this the case, it imposes ncores==1.

Examples

# NOT RUN {
#### 1) Example use
# Set up mixture density
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 <- rmixn(1e4, mu, sigma, w)

# Evaluate density
ds <- dmixn(X, mu, sigma, w = w)
head(ds)

##### 2) More complicated example
# Define mixture
set.seed(5135)
N <- 10000
d <- 2
w <- rep(1, 2) / 2
mu <- matrix(rep(c(0, 0, 2, 3), 2), 2, 2, byrow = TRUE)
sigma <- list(matrix(c(1, 0, 0, 2), 2, 2), matrix(c(1, -0.9, -0.9, 1), 2, 2))

# Simulate random variables
X <- rmixn(N, mu, sigma, w = w, retInd = TRUE)

# Plot mixture density
np <- 100
xvals <- seq(min(X[ , 1]), max(X[ , 1]), length.out = np)
yvals <- seq(min(X[ , 2]), max(X[ , 2]), length.out = np)
theGrid <- expand.grid(xvals, yvals)
theGrid <- as.matrix(theGrid)
dens <- dmixn(theGrid, mu, sigma, w = w)
plot(X, pch = '.', col = attr(X, "index")+1)
contour(x = xvals, y = yvals, z = matrix(dens, np, np),
        levels = c(0.002, 0.01, 0.02, 0.04, 0.08, 0.15 ), add = TRUE, lwd = 2)

# }