projects / morpheus.git / blame
| 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 | |
| 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 | X = rbind(X, newXblock) | |
| 65 | Y = c( Y, vapply(probas, function(p) (rbinom(1,1,p)), 1) ) | |
| 66 | } | |
| 67 | shuffle = sample(n) | |
| 68 | # Returned X should not contain an intercept column (it's an argument of estimation | |
| 69 | # methods) | |
| 70 | list("X"=X[shuffle,-d], "Y"=Y[shuffle], "index"=index[shuffle]) | |
| 71 | } |