mmvstpdfC: C++ implementation of multivariate Normal probability density...

In NPflow: Bayesian Nonparametrics for Automatic Gating of Flow-Cytometry Data

Description

C++ implementation of multivariate Normal probability density function for multiple inputs

Usage

mmvstpdfC(x, xi, psi, sigma, df, Log = TRUE)

Arguments

x	data matrix of dimension p x n, p being the dimension of the data and n the number of data points.
xi	mean vectors matrix of dimension p x K, K being the number of distributions for which the density probability has to be evaluated.
psi	skew parameter vectors matrix of dimension p x K.
sigma	list of length K of variance-covariance matrices, each of dimensions p x p.
df	vector of length K of degree of freedom parameters.
Log	logical flag for returning the log of the probability density function. Defaults is TRUE.

Value

matrix of densities of dimension K x n.

Author(s)

Boris Hejblum

Examples

mmvstpdfC(x = matrix(c(3.399890,-5.936962), ncol=1), xi=matrix(c(0.2528859,-2.4234067)),
psi=matrix(c(11.20536,-12.51052), ncol=1),
sigma=list(matrix(c(0.2134011, -0.0382573, -0.0382573, 0.2660086), ncol=2)),
df=c(7.784106)
)
mvstpdf(x = matrix(c(3.399890,-5.936962), ncol=1), xi=c(0.2528859,-2.4234067),
psi=c(11.20536,-12.51052),
sigma=matrix(c(0.2134011, -0.0382573, -0.0382573, 0.2660086), ncol=2),
df=c(7.784106)
)

#skew-normal limit
mmvsnpdfC(x=matrix(rep(1.96,2), nrow=2, ncol=1),
         xi=matrix(c(0, 0)), psi=matrix(c(1, 1),ncol=1), sigma=list(diag(2))
         )
mvstpdf(x=matrix(rep(1.96,2), nrow=2, ncol=1),
       xi=c(0, 0), psi=c(1, 1), sigma=diag(2),
       df=100000000
       )
mmvstpdfC(x=matrix(rep(1.96,2), nrow=2, ncol=1),
         xi=matrix(c(0, 0)), psi=matrix(c(1, 1),ncol=1), sigma=list(diag(2)),
         df=100000000
         )

#non-skewed limit
mmvtpdfC(x=matrix(rep(1.96,2), nrow=2, ncol=1),
        mean=matrix(c(0, 0)), varcovM=list(diag(2)),
        df=10
        )
mmvstpdfC(x=matrix(rep(1.96,2), nrow=2, ncol=1),
         xi=matrix(c(0, 0)), psi=matrix(c(0, 0),ncol=1), sigma=list(diag(2)),
         df=10
         )

if(require(microbenchmark)){
library(microbenchmark)
microbenchmark(mvstpdf(x=matrix(rep(1.96,2), nrow=2, ncol=1),
                       xi=c(0, 0), psi=c(1, 1),
                       sigma=diag(2), df=10),
               mmvstpdfC(x=matrix(rep(1.96,2), nrow=2, ncol=1),
                         xi=matrix(c(0, 0)), psi=matrix(c(1, 1),ncol=1),
                         sigma=list(diag(2)), df=10),
               times=1000L)
}else{
cat("package 'microbenchmark' not available\n")
}

Example output

Loading required package: Rcpp
Loading required package: truncnorm
          [,1]
[1,] -3.207832
          [,1]
[1,] -3.207832
          [,1]
[1,] -2.986451
          [,1]
[1,] -2.986451
          [,1]
[1,] -2.986451
          [,1]
[1,] -5.258057
          [,1]
[1,] -5.258057
Loading required package: microbenchmark
Unit: microseconds
                                                                                                                                                     expr
                                        mvstpdf(x = matrix(rep(1.96, 2), nrow = 2, ncol = 1), xi = c(0,      0), psi = c(1, 1), sigma = diag(2), df = 10)
 mmvstpdfC(x = matrix(rep(1.96, 2), nrow = 2, ncol = 1), xi = matrix(c(0,      0)), psi = matrix(c(1, 1), ncol = 1), sigma = list(diag(2)),      df = 10)
     min      lq      mean  median       uq      max neval cld
 217.127 243.264 266.47118 256.687 274.6255 3861.909  1000   b
  17.184  21.642  32.82432  31.933  35.5970 1473.457  1000  a 

NPflow documentation built on Feb. 6, 2020, 5:15 p.m.