pkg/R/sampleIO.R

   1 #' Generate sample inputs-outputs
   2 #'
   3 #' Generate input matrix X of size nxd and binary output of size n, where Y is subdivided
   4 #' into K groups of proportions p. Inside one group, the probability law P(Y=1) is
   5 #' described by the corresponding column parameter in the matrix β + intercept b.
   6 #'
   7 #' @name generateSampleIO
   8 #'
   9 #' @param n Number of individuals
  10 #' @param p Vector of K(-1) populations relative proportions (sum (<)= 1)
  11 #' @param β Vectors of model parameters for each population, of size dxK
  12 #' @param b Vector of intercept values (use rep(0,K) for no intercept)
  13 #' @param link Link type; "logit" or "probit"
  14 #'
  15 #' @return A list with
  16 #' \itemize{
  17 #'   \item{X: the input matrix (size nxd)}
  18 #'   \item{Y: the output vector (size n)}
  19 #'   \item{index: the population index (in 1:K) for each row in X}
  20 #' }
  21 #'
  22 #' @examples
  23 #' # K = 3 so we give first two components of p: 0.3 and 0.3 (p[3] = 0.4)
  24 #' io <- generateSampleIO(1000, c(.3,.3),
  25 #'   matrix(c(1,3,-1,1,2,1),ncol=3), c(.5,-1,0), "logit")
  26 #' io$index[1] #number of the group of X[1,] and Y[1] (in 1...K)
  27 #'
  28 #' @export
  29 generateSampleIO <- function(n, p, β, b, link)
  30 {
  31   # Check arguments
  32   tryCatch({n <- as.integer(n)}, error=function(e) stop("Cannot convert n to integer"))
  33   if (length(n) > 1)
  34     warning("n is a vector but should be scalar: only first element used")
  35   if (n <= 0)
  36     stop("n: positive integer")
  37   if (!is.matrix(β) || !is.numeric(β) || any(is.na(β)))
  38     stop("β: real matrix, no NAs")
  39   K <- ncol(β)
  40   if (!is.numeric(p) || length(p)<K-1 || any(is.na(p)) || any(p<0) || sum(p) > 1)
  41     stop("p: positive vector of size >= K-1, no NA, sum(<)=1")
  42   if (length(p) == K-1)
  43     p <- c(p, 1-sum(p))
  44   if (!is.numeric(b) || length(b)!=K || any(is.na(b)))
  45     stop("b: real vector of size K, no NA")
  46
  47   # Random generation of the size of each population in X~Y (unordered)
  48   classes <- rmultinom(1, n, p)
  49
  50   d <- nrow(β)
  51   zero_mean <- rep(0,d)
  52   id_sigma <- diag(rep(1,d))
  53   X <- matrix(nrow=0, ncol=d)
  54   Y <- c()
  55   index <- c()
  56   for (i in 1:K)
  57   {
  58     index <- c(index, rep(i, classes[i]))
  59     newXblock <- MASS::mvrnorm(classes[i], zero_mean, id_sigma)
  60     arg_link <- newXblock %*% β[,i] + b[i]
  61     probas <-
  62       if (link == "logit")
  63       {
  64         e_arg_link = exp(arg_link)
  65         e_arg_link / (1 + e_arg_link)
  66       }
  67       else #"probit"
  68         pnorm(arg_link)
  69     probas[is.nan(probas)] = 1 #overflow of exp(x)
  70     X <- rbind(X, newXblock)
  71     Y <- c( Y, vapply(probas, function(p) (rbinom(1,1,p)), 1) )
  72   }
  73   shuffle <- sample(n)
  74   list("X"=X[shuffle,], "Y"=Y[shuffle], "index"=index[shuffle])
  75 }