simulate data from multinomial logistic regression

simulate-multinomial-LR.R

N <- 1000
J <- 3 # number of predictors (intercept, slp1, slp2)
X <- cbind(1, rnorm(N, 0, 1), rnorm(N, 0, 1)) # N * J
K <- 3 # num of classes
beta <- rbind(c(1, -0.5, 0),
              c(-2, 4, 0),
              c(3, -1.5, 0)) 
# J * K (features * classes) with the last class as the reference
# [cls1_intercept, cls2_intercept, cls3_intercept]
# [cls1_x1,         cls2_x1,        cls3_x1      ]
# [cls1_x2,         cls2_x2,        cls3_x2      ]
z <- X %*% beta
# rowwise softmax transformation
probs <- t(apply(z, 1, function(x) exp(x) / sum(exp(x))))
y <- apply(probs, 1, function(x) rcat(n = 1, prob = x))