tests/testthat/time-test.R

In CubicZebra/MVNBayesian: Bayesian Analysis Framework for MVN (Mixture) Distribution

ipak <- function(pkg) {
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

#libraries
packages <- c("mvtnorm","stats","plyr","Rfast","rgl","commentr")
ipak(packages)

#timecheck
timecheck <- function(code=Sys.sleep(1)){
  start <- Sys.time()
  eval(parse(text=code))
  end <- Sys.time()
  print(end-start)
}

line_comment("check for Ascending_Num()", clipboard = F)
AN <- "x <- c(1,2,2,2,2,2,2,2,3,3,3,1,3,3,3)
x
Ascending_Num(x)"
timecheck(AN)

line_comment("check for MatrixAlternative()", clipboard = F)
MA <- "M <- matrix(1:9,3,3,1)
M
MatrixAlternative(M, 2)"
timecheck(MA)

line_comment("check for MVN_BayesianPosteriori()", clipboard = F)
MVN_BP <- "# Demo using dataset1:
head(dataset1)
BPos <- MVN_BayesianPosteriori(dataset1, c(80,16,3))
BPos$mean
BPos$var

# Singular system caused by insufficient data
eigen(var(dataset1[1:3,]))$values
rcond(var(dataset1[1:3,]))
eigen(var(dataset1[1:6,]))$values
rcond(var(dataset1[1:6,]))

# Singular system caused by improper degree of freedom
K <- cbind(dataset1, dataset1[,3]*(-2)+3)
eigen(var(K[,2:4]))$values
rcond(var(K[,2:4]))"
timecheck(MVN_BP)

line_comment("check for MVN_FConditional()", clipboard = F)
MVN_FC <- "head(dataset1)
BPos <- MVN_BayesianPosteriori(dataset1, c(80,16,3))
BPos # Bayesian Posteriori
result <- MVN_FConditional(BPos, variable = 1, z=c(75, 13, 4))
result$mean
class(result$mean)
result$var
class(result$var)

# compare the following results:
MVN_FConditional(BPos, variable = 2, z=c(75, 13, 4))
MVN_FConditional(BPos, variable = 2, z=c(75, 88, 4))
MVN_FConditional(BPos, variable = 1, z=c(75, 88, 4))"
timecheck(MVN_FC)

line_comment("check for MVN_GibbsSampler()", clipboard = F)
MVN_GS <- "# Get parameters of Bayesian posteriori multivariate normal distribution
BPos <- MVN_BayesianPosteriori(dataset1)
BPos

# Using previous result (BPos) to generate random vectors through Gibbs
# sampling: 7000 observations, start from c(1,1,2), use 0.3 burning rate
BPos_Gibbs <- MVN_GibbsSampler(7000, BPos, initial=c(1,1,2), 0.3)
tail(BPos_Gibbs)

# Check for convergence of Markov chain
BPos$mean
colMeans(BPos_Gibbs)
BPos$var
var(BPos_Gibbs)"
timecheck(MVN_GS)

line_comment("check for MVN_GibbsSampler()", clipboard = F)
############### check for MVN_BayesianIterator() #################
MVN_BI <- "# Bayesian posteriori before iteration using dataset1 as example,
# c(80, 16, 3) as priori mean:
# View 2-norm of covariance matrix of Bayesian posteriori:
BPos_init <- MVN_BayesianPosteriori(dataset1, c(80,16,3))
BPos_init
norm(as.matrix(BPos_init$var), type = \"2\")

# Bayesian posteriori after iteration using c(80,16,3) as priori
# Using 30 last samples generated by GibbsSampler for each step:
BPos_fina1 <- MVN_BayesianIterator(dataset1, c(80,16,3), 5000, 30)
BPos_fina1
norm(as.matrix(BPos_fina1$var), type = \"2\")"
timecheck(MVN_BI)

line_comment("check for MVN_MCMC()", clipboard = F)
MVN_MC <- "# dataset1 has three parameters: fac1, fac2 and fac3:
head(dataset1)

# Get posteriori parameters of dataset1 using prior of c(80,16,3):
BPos <- MVN_BayesianPosteriori(dataset1, pri_mean=c(80,16,3))

# If we want to know when fac1=78, how fac2 responses to fac3, run:
BPos_MCMC <- MVN_MCMC(BPos, steps=8000, pars=c(1), values=c(78), tol=0.3)
MCMC <- BPos_MCMC$MCMCdata
head(MCMC)"
timecheck(MVN_MC)

line_comment("check for MixMVN_BayesianPosteriori()", clipboard = F)
MixMVN_BP <- "# Design matrix should only contain columns of variables
# Export will be a matrix-like data
# Using kmeans (default) clustering algrithm
data_dim <- dataset2[,1:4]
result <- MixMVN_BayesianPosteriori(data=data_dim, species=3)
result

# Get the parameters of the cluster1:
result[1,]

# Get the mixture probability of cluster2:
# (Attention to the difference between
# result[2,1][[1]] and result[2,1])
result[2,1][[1]]

# Get the mean vector of cluster1:
result[1,2][[1]]

# Get the covariance matrix of cluster3:
result[3,3][[1]]"
timecheck(MixMVN_BP)

line_comment("check for MixMVN_GibbsSampler()", clipboard = F)
MixMVN_GS <- "# Use dataset2 for demonstration. Get parameters of Bayesian
# posteriori multivariate normal mixture distribution
head(dataset2)
dataset2_par <- dataset2[,1:4] # only parameter columns are premitted
MixBPos <- MixMVN_BayesianPosteriori(dataset2_par, species=3)
MixBPos

# Generate random vectors using Gibbs sampling:
MixBPos_Gibbs <- MixMVN_GibbsSampler(5000, MixBPos, random_method = \"Gibbs\")
head(MixBPos_Gibbs)

# Compared generation speed of \"Gibbs\" to that of \"Fast\"
MixBPos_Fast <- MixMVN_GibbsSampler(5000, MixBPos, random_method = \"Fast\")
head(MixBPos_Fast)"
timecheck(MixMVN_GS)

line_comment("check for MixMVN_MCMC()", clipboard = F)
MixMVN_MC <- "# dataset2 has 4 parameters: dimen1, dimen2, dimen3 and dimen4:
head(dataset2)
dataset2_dim <- dataset2[,1:4] # extract parametric columns

# Get posteriori parameters of dataset2 using kmeans 3 clustering:
MixBPos <- MixMVN_BayesianPosteriori(dataset2_dim, 3)

# If we want to know when dimen1=1, which clusters are accepted, run:
MixBPos_MCMC <- MixMVN_MCMC(MixBPos, steps=5000, pars=c(1), values=c(1), tol=0.3)
MixBPos_MCMC$AcceptRate
result <- MixBPos_MCMC$MCMCdata
head(result)

# count accepted samples by clustering:
count(result[which(result[,7]==1),5])"
timecheck(MixMVN_MC)

rm(list = ls())