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