update
[morpheus.git] / pkg / R / sampleIO.R
... / ...
CommitLineData
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
29generateSampleIO <- 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}