Hello Stan users and exports,

I have a problem as the tittle describes.

My model is a little complex so I paste the parts that cause error.

```
data{
int<lower=0> y[3,18];
real log_max; // log(1000)
//etc...
}
parameters{
vector<lower = 0>[2] sigma;
vector[2] U;
matrix[2,2] B;
matrix<upper=log_max>[2,18]log_y_2;
//etc...
}
transformed parameters{
simplex[3] theta[18];
matrix<upper=log_max>[3,18]log_y;
for(t in 1:18){
log_y[1,t]=log_y_2[1,t];
log_y[2,t]=log_y_2[2,t];
if ((exp(log_max)-exp(log_y_2[1,t])-exp(log_y_2[2,t]))<25)
log_y[3,t]=log(25);
else
log_y[3,t]=log(exp(log_max)-exp(log_y_2[1,t])-exp(log_y_2[2,t]));
}
for (t in 1:18){
theta[t]=softmax(exp(log_y[,t]));
}
}
model{
for (i in 1:18) {
y[,i] ~ multinomial(theta[i,]);
}
for(t in 2:18){
target+= multi_normal_lpdf(log_y_2[,t] |U+B*log_y_2[,t-1],diag_matrix(sigma));
//etc...
}
```

After many attempts, I found that log_y_3[3,t] in transformed parameters is the problem. If I write

```
if ((exp(log_max)-exp(log_y_2[1,t])-exp(log_y_2[2,t]))<25)
log_y[3,t]=log(25);
else
log_y[3,t]=log(exp(log_max)-exp(log_y_2[1,t])-exp(log_y_2[2,t]));
```

there would be a log(0) error.

However, if I write

```
if ((exp(log_max)-exp(log_y_2[1,t])-exp(log_y_2[2,t]))<25)
log_y[3,t]=log(25);
else
log_y[3,t]=log(50);
```

It will work.

But if change 50 to 1000 or 500, it will report an error once again.

Actually, I need **log_y[3,t]=log(exp(log_max)-exp(log_y_2[1,t])-exp(log_y_2[2,t]));**. Because log_y is a matrix that elements in each column follow multinomial distribution. How can I successfully achieve this goal?

Thank you in advance