# Multinomial regression with restriction

Dear Stan Users,

I am new in stan and I would like to fit a multinomial model with a restriction on the probabilities. I believe it is different from a softmax model, but I am not really sure about it.
Below, you can see my stan model and R code. However, when compiling in R a got some error messages. I looked for help on internet, but no success so far. Please, let me know if some of you already work with a similar model and how to code it in stan.

Thanks.

``````data {
int K;
int N;
int D;
int y[N, K];
matrix[N, D] X;
}
parameters {
real beta1;
real beta2;
}
transformed parameters {
real linpred1[N];
real linpred2[N];
matrix[N, D] prob;
for(i in 1:N) {
linpred1[i] = exp(X[i,1]*beta1 + X[i,2]*beta1 + X[i,3]*beta1);
linpred2[i] = exp(X[i,1]*beta2 + X[i,2]*beta2 + X[i,3]*beta2);
prob[i,2] = linpred1[i]/(1+linpred1[i] + linpred2[i]);
prob[i,3] = linpred2[i]/(1+linpred1[i] + linpred2[i]);
prob[i,1] = 1 - prob[i,2] - prob[i,3];
}
}
model {
for (i in 1:N) {
y[i,] ~ multinomial(to_row_vector(prob[i,]) );
}
}
``````
``````set.seed(123)

# Generating the probabilities ---------------------------------------
x1 <- seq(-1, 1, l = 100)
x2 <- rbinom(100, prob = 0.5, size = 1)
X <- model.matrix(~ x1 + x2)
beta1 <- c(0, 0.8, -1)
beta2 <- c(0.5, -0.8, -1)

eta1 <- X%*%beta1
eta2 <- X%*%beta2

pi1 <- exp(eta1)/(1+ exp(eta1) + exp(eta2))
pi2 <- exp(eta2)/(1+ exp(eta1) + exp(eta2))
pi0 <- 1 - pi1 - pi2

# Simulating data set ------------------------------------------------
Y <- matrix(NA, ncol = 3, nrow = 100)
for(i in 1:100) {
Y[i,] <- rmultinom(1, size = 1, prob = cbind(pi0, pi1, pi2)[i,])
}

library(rstan)
data <- list("y" = Y, "K" = 3, "D" = 3, "N" = 100, "X" = X,
beta1 = c(0.1, 0.6, 0), beta2 = c(0.1, 0.8, 0))

fit1 <- stan(
file = "mult1.stan",  # Stan program
data = data,    # named list of data
chains = 4,             # number of Markov chains
warmup = 1000,          # number of warmup iterations per chain
iter = 2000,            # total number of iterations per chain
cores = 2,              # number of cores (could use one per chain)
refresh = 0             # no progress shown
)
``````