Movatterモバイル変換

Fast density computation for mixture of multivariate normal distributions.

Description

Fast density computation for mixture of multivariate normal distributions.

Usage

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

Arguments

Details

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

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 ofX).

Author(s)

Matteo Fasiolo <matteo.fasiolo@gmail.com>.

Examples

#### 1) Example use# Set up mixture densitymu <- matrix(c(1, 2, 10, 20), 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)# SimulateX <- rmixn(1e4, mu, sigma, w)# Evaluate densityds <- dmixn(X, mu, sigma, w = w)head(ds)##### 2) More complicated example# Define mixtureset.seed(5135)N <- 10000d <- 2w <- rep(1, 2) / 2mu <- matrix(c(0, 0, 2, 3), 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 variablesX <- rmixn(N, mu, sigma, w = w, retInd = TRUE)# Plot mixture densitynp <- 100xvals <- 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)

Fast density computation for mixture of multivariate Student's t distributions.

Description

Fast density computation for mixture of multivariate Student's t distributions.

Usage

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

Arguments

Details

There are many candidates for the multivariate generalization of Student's t-distribution, here we usethe parametrization described herehttps://en.wikipedia.org/wiki/Multivariate_t-distribution. NB: at the moment the parallelization does not work properly on Solaris OS whenncores>1. Hence,dmixt() checks if the OS is Solaris and, if this the case, it imposesncores==1.

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 ofX).

Author(s)

Matteo Fasiolo <matteo.fasiolo@gmail.com>.

Examples

#### 1) Example use# Set up mixture densitydf <- 6mu <- matrix(c(1, 2, 10, 20), 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)# SimulateX <- rmixt(1e4, mu, sigma, df, w)# Evaluate densityds <- dmixt(X, mu, sigma, w = w, df = df)head(ds)##### 2) More complicated example# Define mixtureset.seed(5135)N <- 10000d <- 2df = 10w <- rep(1, 2) / 2mu <- matrix(c(0, 0, 2, 3), 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 variablesX <- rmixt(N, mu, sigma, w = w, df = df, retInd = TRUE)# Plot mixture densitynp <- 100xvals <- 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 <- dmixt(theGrid, mu, sigma, w = w, df = df)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)

Fast computation of the multivariate normal density.

Description

Fast computation of the multivariate normal density.

Usage

dmvn(X, mu, sigma, log = FALSE, ncores = 1, isChol = FALSE)

Arguments

Value

A vector of length n where the i-the entry contains the pdf of the i-th random vector.

Author(s)

Matteo Fasiolo <matteo.fasiolo@gmail.com>

Examples

N <- 100d <- 5mu <- 1:dX <- t(t(matrix(rnorm(N*d), N, d)) + mu)tmp <- matrix(rnorm(d^2), d, d)mcov <- tcrossprod(tmp, tmp)  + diag(0.5, d)myChol <- chol(mcov)head(dmvn(X, mu, mcov), 10)head(dmvn(X, mu, myChol, isChol = TRUE), 10)## Not run: # Performance comparison: microbenchmark does not work on all# platforms, hence we need to check whether it is installedif( "microbenchmark" %in% rownames(installed.packages()) ){library(mvtnorm)library(microbenchmark)a <- cbind(      dmvn(X, mu, mcov),      dmvn(X, mu, myChol, isChol = TRUE),      dmvnorm(X, mu, mcov))      # Check if we get the same output as dmvnorm()a[ , 1] / a[, 3]a[ , 2] / a[, 3]microbenchmark(dmvn(X, mu, myChol, isChol = TRUE),                dmvn(X, mu, mcov),                dmvnorm(X, mu, mcov))               detach("package:mvtnorm", unload=TRUE)}## End(Not run)

Fast computation of the multivariate Student's t density.

Description

Fast computation of the multivariate Student's t density.

Usage

dmvt(X, mu, sigma, df, log = FALSE, ncores = 1, isChol = FALSE)

Arguments

Details

There are many candidates for the multivariate generalization of Student's t-distribution, here we usethe parametrization described herehttps://en.wikipedia.org/wiki/Multivariate_t-distribution. NB: at the moment the parallelization does not work properly on Solaris OS whenncores>1. Hence,dmvt() checks if the OS is Solaris and, if this the case, it imposesncores==1.

Value

A vector of length n where the i-the entry contains the pdf of the i-th random vector.

Author(s)

Matteo Fasiolo <matteo.fasiolo@gmail.com>

Examples

N <- 100d <- 5mu <- 1:ddf <- 4X <- t(t(matrix(rnorm(N*d), N, d)) + mu)tmp <- matrix(rnorm(d^2), d, d)mcov <- tcrossprod(tmp, tmp)  + diag(0.5, d)myChol <- chol(mcov)head(dmvt(X, mu, mcov, df = df), 10)head(dmvt(X, mu, myChol, df = df, isChol = TRUE), 10)

Fast computation of squared mahalanobis distance between all rows ofX and the vectormu with respect to sigma.

Description

Fast computation of squared mahalanobis distance between all rows ofX and the vectormu with respect to sigma.

Usage

maha(X, mu, sigma, ncores = 1, isChol = FALSE)

Arguments

Value

a vector of length n where the i-the entry contains the square mahalanobis distance i-th random vector.

Author(s)

Matteo Fasiolo <matteo.fasiolo@gmail.com>

Examples

N <- 100d <- 5mu <- 1:dX <- t(t(matrix(rnorm(N*d), N, d)) + mu)tmp <- matrix(rnorm(d^2), d, d)mcov <- tcrossprod(tmp, tmp)myChol <- chol(mcov)rbind(head(maha(X, mu, mcov), 10),      head(maha(X, mu, myChol, isChol = TRUE), 10),      head(mahalanobis(X, mu, mcov), 10))## Not run: # Performance comparison: microbenchmark does not work on all# platforms, hence we need to check whether it is installedif( "microbenchmark" %in% rownames(installed.packages()) ){library(microbenchmark)a <- cbind(  maha(X, mu, mcov),  maha(X, mu, myChol, isChol = TRUE),  mahalanobis(X, mu, mcov))  # Same output as mahalanobisa[ , 1] / a[, 3]a[ , 2] / a[, 3]microbenchmark(maha(X, mu, mcov),               maha(X, mu, myChol, isChol = TRUE),               mahalanobis(X, mu, mcov))}## End(Not run)

Mean-shift mode seeking algorithm

Description

Given a sample from a d-dimensional distribution, an initialization point and a bandwidththe algorithm finds the nearest mode of the corresponding Gaussian kernel density.

Usage

ms(X, init, H, tol = 1e-06, ncores = 1, store = FALSE)

Arguments

Value

A list whereestim is a d-dimensional vector containing the last position of the algorithm, whiletraj is a matrix with d-colums representing the trajectory of the algorithm along each dimension. Ifstore == FALSE the whole trajectoryis not stored andtraj = NULL.

Author(s)

Matteo Fasiolo <matteo.fasiolo@gmail.com>.

Examples

set.seed(434)# Simulating multivariate normal dataN <- 1000mu <- c(1, 2)sigma <- matrix(c(1, 0.5, 0.5, 1), 2, 2)X <- rmvn(N, mu = mu, sigma = sigma)# Plotting the true density functionsteps <- 100range1 <- seq(min(X[ , 1]), max(X[ , 1]), length.out = steps)range2 <- seq(min(X[ , 2]), max(X[ , 2]), length.out = steps)grid <- expand.grid(range1, range2)vals <- dmvn(as.matrix(grid), mu, sigma)contour(z = matrix(vals, steps, steps),  x = range1, y = range2, xlab = "X1", ylab = "X2")points(X[ , 1], X[ , 2], pch = '.') # Estimating the mode from "nrep" starting pointsnrep <- 10index <- sample(1:N, nrep)for(ii in 1:nrep) {  start <- X[index[ii], ]  out <- ms(X, init = start, H = 0.1 * sigma, store = TRUE)  lines(out$traj[ , 1], out$traj[ , 2], col = 2, lwd = 2)   points(out$final[1], out$final[2], col = 4, pch = 3, lwd = 3) # Estimated mode (blue)  points(start[1], start[2], col = 2, pch = 3, lwd = 3)         # ii-th starting value }

Fast simulation of r.v.s from a mixture of multivariate normal densities

Description

Fast simulation of r.v.s from a mixture of multivariate normal densities

Usage

rmixn(  n,  mu,  sigma,  w,  ncores = 1,  isChol = FALSE,  retInd = FALSE,  A = NULL,  kpnames = FALSE)

Arguments

Details

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 thisRNG 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 theset.seed() function can be used to obtain reproducible results. Notice though that changingncores causes most of the generated numbersto be different even if R's seed is the same (see example below). NB: at the moment the RNG does not workproperly on Solaris OS whenncores>1. Hence,rmixn() checks if the OS is Solaris and, if this the case, it imposesncores==1.

Value

IfA==NULL (default) the output is an (n x d) matrix where the i-th row is the i-th simulated vector.IfA!=NULL then the random vector are store inA, which is provided by the user, and the functionreturnsNULL. Notice that ifretInd==TRUE an attribute called "index" will be added to A.This is a vector specifying to which mixture components each random vector belongs.

Author(s)

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 componentsmu <- matrix(c(1, 2, 10, 20), 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)# SimulateX <- rmixn(1e4, mu, sigma, w, retInd = TRUE)plot(X, pch = '.', col = attr(X, "index"))# Simulate with fixed seedset.seed(414)rmixn(4, mu, sigma, w)set.seed(414)rmixn(4, mu, sigma, w)set.seed(414)  rmixn(4, mu, sigma, 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)rmixn(4, mu, sigma, w, ncores = 2, A = A) # This returns NULL ...A                                         # ... but the result is here

Fast simulation of r.v.s from a mixture of multivariate Student's t densities

Description

Fast simulation of r.v.s from a mixture of multivariate Student's t densities

Usage

rmixt(  n,  mu,  sigma,  df,  w,  ncores = 1,  isChol = FALSE,  retInd = FALSE,  A = NULL,  kpnames = FALSE)

Arguments

Details

There are many candidates for the multivariate generalization of Student's t-distribution, here we usethe parametrization described herehttps://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 thisRNG 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 theset.seed() function can be used to obtain reproducible results. Notice though that changingncores causes most of the generated numbersto 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 whenncores>1. Hence,rmixt() checks if the OS is Solaris and, if this the case, it imposesncores==1

Value

IfA==NULL (default) the output is an (n x d) matrix where the i-th row is the i-th simulated vector.IfA!=NULL then the random vector are store inA, which is provided by the user, and the functionreturnsNULL. Notice that ifretInd==TRUE an attribute called "index" will be added to A.This is a vector specifying to which mixture components each random vector belongs.

Author(s)

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 componentsdf <- 6mu <- matrix(c(1, 2, 10, 20), 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)# SimulateX <- rmixt(1e4, mu, sigma, df, w, retInd = TRUE)plot(X, pch = '.', col = attr(X, "index"))# Simulate with fixed seedset.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

Fast simulation of multivariate normal random variables

Description

Fast simulation of multivariate normal random variables

Usage

rmvn(n, mu, sigma, ncores = 1, isChol = FALSE, A = NULL, kpnames = FALSE)

Arguments

Details

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 thisRNG 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 theset.seed() function can be used to obtain reproducible results. Notice though that changingncores causes most of the generated numbersto be different even if R's seed is the same (see example below). NB: at the moment the RNG does not workproperly on Solaris OS whenncores>1. Hence,rmvn() checks if the OS is Solaris and, if this the case, it imposesncores==1.

Value

IfA==NULL (default) the output is an (n x d) matrix where the i-th row is the i-th simulated vector.IfA!=NULL then the random vector are store inA, which is provided by the user, and the functionreturnsNULL.

Author(s)

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

d <- 5mu <- 1:d# Creating covariance matrixtmp <- matrix(rnorm(d^2), d, d)mcov <- tcrossprod(tmp, tmp)set.seed(414)rmvn(4, 1:d, mcov)set.seed(414)rmvn(4, 1:d, mcov)set.seed(414)  rmvn(4, 1:d, mcov, 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, d)class(A) <- "numeric" # This is important. We need the elements of A to be of class "numeric". set.seed(414)rmvn(4, 1:d, mcov, ncores = 2, A = A) # This returns NULL ...A                                     # ... but the result is here

Fast simulation of multivariate Student's t random variables

Description

Fast simulation of multivariate Student's t random variables

Usage

rmvt(n, mu, sigma, df, ncores = 1, isChol = FALSE, A = NULL, kpnames = FALSE)

Arguments

Details

There are in fact many candidates for the multivariate generalization of Student's t-distribution, here we usethe parametrization described herehttps://en.wikipedia.org/wiki/Multivariate_t-distribution.

Notice thatrmvt() 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 thisRNG 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 theset.seed() function can be used to obtain reproducible results. Notice though that changingncores causes most of the generated numbersto be different even if R's seed is the same (see example below). NB: at the moment the RNG does not workproperly on Solaris OS whenncores>1. Hence,rmvt() checks if the OS is Solaris and, if this the case, it imposesncores==1.

Value

IfA==NULL (default) the output is an (n x d) matrix where the i-th row is the i-th simulated vector.IfA!=NULL then the random vector are store inA, which is provided by the user, and the functionreturnsNULL.

Author(s)

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

d <- 5mu <- 1:ddf <- 4# Creating covariance matrixtmp <- matrix(rnorm(d^2), d, d)mcov <- tcrossprod(tmp, tmp) + diag(0.5, d)set.seed(414)rmvt(4, 1:d, mcov, df = df)set.seed(414)rmvt(4, 1:d, mcov, df = df)set.seed(414)  rmvt(4, 1:d, mcov, df = df, ncores = 2) # These will not match the r.v. generated on a single core.###### Here we create the matrix that will hold the simulated random variables upfront.A <- matrix(NA, 4, d)class(A) <- "numeric" # This is important. We need the elements of A to be of class "numeric". set.seed(414)rmvt(4, 1:d, mcov, df = df, ncores = 2, A = A) # This returns NULL ...A                                     # ... but the result is here