Hi all,

I’m essentially trying to sample from a Bayesian multinomial logit with a few complexities: choices have characteristics and choices might not be available in all time periods (no time series here). Choice probabilities look sth like `s_jt = softmax(d + X * beta)`

adjusted for availability. X is the feature matrix of size J products * L features, which includes a constant so I use priors:

- for variance of d:
`Omega = LKJ( J, 2.0)`

, thus`d ~ MvNormal(0.0, Omega)`

- for beta:
`beta ~ MvNormal(ones(L), diagm(10 * ones(L)))`

I then loop over time periods to compute the choice probabilities with sth as such:

```
for t in 1:T
observed_qjt[:, t] ~ Multinomial(mkt_size, sjt)
end
```

The code seems to work in the sense that I don’t get any errors and the sampling procedure begins, but it stays stuck at 0% even when I reduce the number of time periods. To be precise about dimensions of the problem, I have J = 45 choices, T = 1024 time periods and L = 4 features.

I have sth that’s close to a MWE if you need but the matrices are quite big so not sure what’s the best way to share them.

```
mAV = BitMatrix(vcat(data[:mAV], ones(1, size(data[:mAV], 2))))[:, mkts]
mX = Matrix{Float64}(data[:mX])
mQ = Matrix{Int64}(round.(data[:mQ]))[:, mkts]
# Priors
J = size(mX, 1)
T = size(mQ, 2)
β₀ = ones(size(mX, 2))
Vᵦ = diagm(10 * β₀)
# Settings for the Hamiltonian Monte Carlo HMC sampler
iterations = 1_000 #Number of samples
ϵ = 0.1
τ = 5
@model function estimateBayesSharesDiv_logit(mAV, mX, mQ, J, T, β₀, Vᵦ)
# Instantiate intermediate objects
δ = Matrix{Real}(undef, J, T)
# Priors
Ω_prod ~ LKJ(J, 2.0) # on δⱼ
β ~ MvNormal(β₀, Vᵦ) # L chars
# Define parameters
δ .~ MvNormal(zeros(J), Symmetric(Ω_prod))
# Market specific probabilities
for t in 1:T
# Extract availability
av = mAV[:, t]
# Instantiate data holders
u_jt = zeros(eltype(mX[1][1] * β), J+1)
s_jt = zeros(eltype(mX[1][1] * β), J+1)
# Compute utility
u_jt[1:J] = mX * β + δ[:, t]
u_jt[end] = 0.0
# Compute choice probabilities for available prods
s_jt[av] = softmax(u_jt[av])
# Data ~ Model
mkt_size = Int(round(sum(mQ[:, t])))
mQ[:, t] ~ Multinomial(mkt_size, s_jt)
end
return mQ
end
# Sampling
# chn = sample(estimateBayesSharesDiv_logit(mAV, mX, mQ, J, T, β₀, Vᵦ), HMC(ϵ, τ), iterations)
chn = sample(estimateBayesSharesDiv_logit(mAV, mX, mQ, J, T, β₀, Vᵦ), NUTS(0.65), iterations)
```