experiments/Table01-10/Table3.R

In hkasahar/normalregMix: normalregMix: An R Package for Normal Mixture Regression Models

## 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)
}

musets <- list(c(-4.5,-1.5, 1.5, 4.5), c(-6,-2, 2, 6),
               c(-4.0,-1.25, 1.25, 4.0), c(-3.5,-1.5, 1.5, 5.0))
sigmasets <- list(c(1, 1, 1, 1), c(0.6, 1.2, 0.6, 1.2),
                  c(0.6, 0.8, 1.0, 1.2), c(0.8, 0.8, 0.8, 1.2))
alphasets <- list(rep(1/4, 4))
nset <- c(200, 400)
table03.pairs <- GeneratePhiDataPairs(
  CreatePhiSet(musets, sigmasets, alphasets, nset))

## 3. Calling Tables
# Table 03
print("Table 03:")
GenerateTable0203(table03.pairs)

## ====== END EXPERIMENT ======
# Rmpi termination
mpi.close.Rslaves()