## PARAMS
RMPI <- FALSE
REPLICATION <- 5
## Generates EM test result according to the dimension of X
PerformEMtest <- function (sample, q, m = 1, z = NULL, parallel) {
library(doParallel) # workers might need information
library(normalregMix) # workers might need information
n <- as.integer(length(sample)/(q+1))
y <- sample[1:n] # first n elements represents y data
if (q <= 0)
return (normalmixMEMtest(y, m = m, z = z, crit.method = "asy",
parallel = parallel))
# the other part consists of n chuck of q-length x data
x <- matrix(sample[(n+1):length(sample)], nrow = n, byrow = TRUE)
return (regmixMEMtest(y, x, m = m, z = z, crit.method = "asy",
parallel = parallel))
}
## Returns frequency that the null H0: m=1 is rejected
# out of replications of given an and data that consists of columns of samples
PerformEMtests <- function (data, crit = 0.05, q = 0, m = 1,
parallel = 1) {
if (RMPI)
{
# need to transform data (matrix) to a list first; each element is a column (y x_1' x_2' ... x_n')'
ldata <- lapply(seq_len(ncol(data)), function(i) data[,i])
out <- mpi.applyLB(ldata, PerformEMtest, q = q, m = m, z = NULL,
parallel = 0)
}
else
out <- apply(data, 2, PerformEMtest, q = q, m = m, z = NULL,
parallel = parallel)
pvals <- sapply(out, "[[", "pvals")
return (list(K1 = mean(pvals[1,] < crit),
K2 = mean(pvals[2,] < crit),
K3 = mean(pvals[3,] < crit)))
}
# Generate a column that represents a sample using phi given.
GenerateSample <- function(phi, dim.x){
n <- phi$n
alpha <- phi$alphaset
mus <- phi$muset
sigmas <- phi$sigmaset
m <- length(alpha)
alpha_cumsum <- cumsum(alpha)
betaset <- matrix(rep(0,m), ncol = m)
sample.x <- matrix(rep(0,n), ncol = 1)
sample.y <- vector(mode = "double", n)
if (!is.null(phi$betaset))
{
if (dim.x == 0) # if dim.x is given zero, it is assumed dim.x = 1.
dim.x = 1
betaset <- t(replicate(dim.x, unlist(phi$betaset)))
sample.x <- matrix(rnorm(n*dim.x), nrow = n) # each row is one observation
}
for (i in 1:n)
{
prob <- runif(1)
index <- 1
for (j in 2:m)
if (prob > alpha_cumsum[j-1] && prob <= alpha_cumsum[j])
index = j
sample.y[i] <- rnorm(1, mean = (mus[index] + betaset[,index]*sample.x[i,]),
sd = sigmas[index])
}
if (dim.x > 0)
return (rbind(sample.y, (t(sample.x))))
return(sample.y)
}
## Generate a pair of phi and data, where data is generated by replication.
GeneratePhiDataPair <- function(phi, replication, dim.x) {
phi <- as.list(phi) # make it atomic.
# data is an (n replication) matrix whose column represents a sample of size n,
data <- do.call(cbind,
replicate(replication,
GenerateSample(phi = phi, dim.x = dim.x), simplify = FALSE))
return (list(phi = phi, data = data))
}
## Create data given phiset and replication
# returns a list of (phi = phi, data = data)
GeneratePhiDataPairs <- function(phiset, replication = 10, dim.x = 0) {
apply(phiset, 1, GeneratePhiDataPair, replication = replication, dim.x = dim.x)
}
## Returns a list of phis given musets, sigmasets, alphasets, nset, betasets
CreatePhiSet <- function(musets, sigmasets, alphasets, nset, betasets = NULL) {
if (is.null(betasets))
return (expand.grid(muset = musets, sigmaset = sigmasets, alphaset = alphasets,
n = nset))
expand.grid(muset = musets, sigmaset = sigmasets, alphaset = alphasets,
n = nset, betaset = betasets)
}
### EXPERIMENT
## Rmpi setup
print("collecting workers..")
mpi.spawn.Rslaves()
mpi.setup.rngstream()
mpi.bcast.Robj2slave(PerformEMtest, all=TRUE)
print("workers loaded.")
## ====== BEGIN EXPERIMENT ======
## 1. Initialization & data generation
set.seed(123456)
# Table 02
musets <- list(c(-2.0, 0, 2.0), c(-1.0, 0, 3.0))
sigmasets <- list(c(0.6, 1.2, 0.6))
alphasets <- list(c(1/3, 1/3, 1/3), c(0.4, 0.2, 0.4))
phi1 <- CreatePhiSet(musets, sigmasets, alphasets, nset = 200)
phi2 <- CreatePhiSet(musets, sigmasets, alphasets, nset = 400)
table02.pairs <- GeneratePhiDataPairs(rbind(phi1, phi2),
replication = REPLICATION)
## 2. Tables
GenerateTable0203 <- function(pairs)
{
rows <- list()
for (i in 1:length(pairs)) {
phi <- pairs[[i]]$phi
data <- pairs[[i]]$data
m <- length(phi$alphaset)
powers <- PerformEMtests(data = data, m = (m-1), crit = 0.05)
rows[[i]] <- list(phi$muset, phi$sigmaset, phi$alphaset, phi$n,
powers$K1, powers$K2, powers$K3)
df <- data.frame(matrix(unlist(rows), ncol = length(unlist(rows[[1]])), byrow=T))
colnames(df) <- c(paste("mu",1:m,sep=''), paste("sigma",1:m,sep=''), paste("alpha",1:m,sep=''),
"n", "pow.K1", "pow.K2", "pow.K3")
print(df) # save every time
}
print(df)
}
## 3. Calling Tables
# Table 02
print("Table 02:")
GenerateTable0203(table02.pairs)
## ====== END EXPERIMENT ======
# Rmpi termination
mpi.close.Rslaves()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.