| 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 |
| 8 | #' @param p Vector of K-1 populations relative proportions (sum <= 1) |
| 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 | #' |
| 20 | #' @export |
| 21 | generateSampleIO = function(n, p, β, b, link) |
| 22 | { |
| 23 | # Check arguments |
| 24 | tryCatch({n = as.integer(n)}, error=function(e) stop("Cannot convert n to integer")) |
| 25 | if (length(n) > 1) |
| 26 | warning("n is a vector but should be scalar: only first element used") |
| 27 | if (n <= 0) |
| 28 | stop("n: positive integer") |
| 29 | if (!is.matrix(β) || !is.numeric(β) || any(is.na(β))) |
| 30 | stop("β: real matrix, no NAs") |
| 31 | K = ncol(β) |
| 32 | if (!is.numeric(p) || length(p)!=K-1 || any(is.na(p)) || any(p<0) || sum(p) > 1) |
| 33 | stop("p: positive vector of size K-1, no NA, sum<=1") |
| 34 | p <- c(p, 1-sum(p)) |
| 35 | if (!is.numeric(b) || length(b)!=K || any(is.na(b))) |
| 36 | stop("b: real vector of size K, no NA") |
| 37 | |
| 38 | #random generation of the size of each population in X~Y (unordered) |
| 39 | classes = rmultinom(1, n, p) |
| 40 | |
| 41 | d = nrow(β) |
| 42 | zero_mean = rep(0,d) |
| 43 | id_sigma = diag(rep(1,d)) |
| 44 | # Always consider an intercept (use b=0 for none) |
| 45 | d = d + 1 |
| 46 | β = rbind(β, b) |
| 47 | X = matrix(nrow=0, ncol=d) |
| 48 | Y = c() |
| 49 | index = c() |
| 50 | for (i in 1:ncol(β)) |
| 51 | { |
| 52 | index = c(index, rep(i, classes[i])) |
| 53 | newXblock = cbind( MASS::mvrnorm(classes[i], zero_mean, id_sigma), 1 ) |
| 54 | arg_link = newXblock %*% β[,i] #β |
| 55 | probas = |
| 56 | if (link == "logit") |
| 57 | { |
| 58 | e_arg_link = exp(arg_link) |
| 59 | e_arg_link / (1 + e_arg_link) |
| 60 | } |
| 61 | else #"probit" |
| 62 | pnorm(arg_link) |
| 63 | probas[is.nan(probas)] = 1 #overflow of exp(x) |
| 64 | #probas = rowSums(p * probas) |
| 65 | X = rbind(X, newXblock) |
| 66 | #Y = c( Y, vapply(probas, function(p) (ifelse(p >= .5, 1, 0)), 1) ) |
| 67 | Y = c( Y, vapply(probas, function(p) (rbinom(1,1,p)), 1) ) |
| 68 | } |
| 69 | shuffle = sample(n) |
| 70 | # Returned X should not contain an intercept column (it's an argument of estimation |
| 71 | # methods) |
| 72 | list("X"=X[shuffle,-d], "Y"=Y[shuffle], "index"=index[shuffle]) |
| 73 | } |