I am trying to fit a multinomial probit model via data augmentation with Stan where the predictors are fixed effects for each dimension, and one alternative-specific variable that varies in the cross-section. In practical terms, I want to model the demand for products via their latent attributes (fixed effects) and the price that is unique to every customer.

For simplicity, I’ve fixed the covariance matrix at an identity, such that only the coefficients (ie preferences) are to be estimated. With simulated data, the model does however not recover the parameters, and it is unclear to me why. The comments by @bgoodri and @betanalpha in [this](// Scaling of Log-likelihood value in Hamiltonian Monte Carlo - #5 by Subodh) post I am aware of, but (perhaps naively) I assumed that fixing the covariance matrix at an identity could perhaps alleviate the issues that other people have had.

What is particularly puzzling to me is that the E-BFMI is generally between 0.3 to 0.5, and the other diagnostics (from my understanding of them) are not particularly irregular either.

If anyone has a deep understanding of this, could you point me to potential errors in my model specification or explain why the model cannot be fitted properly?

Here’s my code:

```
functions {
real row_max_except(row_vector a, int d) {
int D = cols(a);
real max_val;
if (d == 1) {
max_val = max(a[2 : D]);
} else if (d == D) {
max_val = max(a[1 : (D - 1)]);
} else {
max_val = max(a[1 : (d - 1)]);
max_val = fmax(max_val, max(a[(d + 1) : D]));
}
return fmax(max_val, 0);
}
real partial_sum_lpdf(array[] matrix comb, int start, int end, vector beta,
matrix Sigma, int D) {
array[end - start + 1] vector[D] xb;
for (i in 1 : (end - start + 1)) {
xb[i] = comb[i, : , 2 : (D + 2)] * beta;
}
return multi_normal_lupdf(comb[ : , : , 1] | xb, Sigma);
}
}
data {
int<lower=1> D; // Number of dimensions
int<lower=1> N; // Number of observations
array[N] int<lower=0, upper=D + 1> y; // Outcome category
matrix[N * D, D + 1] x; // Design matrix: Stacked identities of size D and a column vector of a price variable
}
parameters {
vector[D + 1] beta; // Regression coefficients
matrix[N, D] z_raw; // Unconstrained utilities
}
transformed parameters {
matrix[N, D] z; // Constrained utilities
profile("utilconst") {
for (n in 1 : N) {
for (d in 1 : D) {
if (y[n] == d) {
z[n, d] = exp(z_raw[n, d]) + row_max_except(row(z_raw, n), d);
} else {
z[n, d] = -exp(z_raw[n, d]) + row_max_except(row(z_raw, n), d);
}
}
}
}
array[N] matrix[D, D + 2] comb;
profile("combine") {
for (n in 1 : N) {
comb[n] = append_col(z[n, : ]', x[((n - 1) * D + 1) : (n * D), : ]);
}
}
}
model {
profile("params") {
beta ~ normal(0, 1);
}
profile("choice") {
matrix[D, D] Sigma = identity_matrix(D);
target += reduce_sum(partial_sum_lupdf, comb, 1, beta, Sigma, D);
}
profile("jacadj") {
for (n in 1 : N) {
for (d in 1 : D) {
target += z_raw[n, d]; // Jacobian is exp(z_raw) in both cases, adjustment is the log Jacobian
}
}
}
}
```