I am trying to use Turing to estimate an ordinal multinomial model with partially fixed thresholds/cutpoints, but I’m wondering if I’m doing something wrong: metropolis-hastings sampling works great, but the autodiff-based samplers (NUTS/HMC) can’t seem to do a single iteration, even running overnight (plus indicate errors in estimating starting values for the sampler).

Any Turing users have any ideas what’s going wrong here/what I’m doing wrong?

I’m wondering if it’s a kind of type instability (perhaps from combining the fixed cutpoints with the estimated ones into a single vector), but (a) I was having a lot of trouble finding correct syntax to check that in current versions of Turing, and (b) even if that’s the problem, I still wouldn’t know how to solve it.

EDIT: Just after I posted, I found the correct code for checking type instability. The code involving the multinomial distribution and the ps vector show up in red.

Here’s a MWE:

```
using Turing, Distributions
responses = [1, 5, 23, 21, 27] # Number of people choosing responses 1-5
total_n = 77
@model function ordered_multinomial(responses, counts, threshold_lower, threshold_upper)
# model mean and std of the latent normal distribution
mu ~ Normal(0, 2)
sigma ~ Exponential(2)
# Model cutpoints/thresholds. There are K-1 thresholds
# (where K is the number of response options).
# Here, two thresholds are fixed (to allow estimating the
# latent normal location and scale). To estimate the two
# "free" thresholds, I partition the space between the upper
# and lower thresholds into three parts, using a Dirichlet prior
threshold_partitions ~ Dirichlet(3, 1)
# Make a vector of all thresholds/cutpoints, combining known
# (fixed) values and estimated values
threshold_points = cumsum(threshold_partitions) * (threshold_upper - threshold_lower) .+ threshold_lower
c = vcat(-Inf, threshold_lower, threshold_points, Inf)
# Using the latent normal distribution and set of cutpoints, get probabilities of
# people choosing a given response option (1-5)
ps = [cdf(Normal(0, 1), (i[2] - mu) / sigma) - cdf(Normal(0, 1), (i[1] - mu) / sigma ) for i in zip(c[1:(end-1)], c[2:end])]
# Model likelihood of response pattern given set of response
# probabilities and a total sample size
responses ~ Multinomial(counts, ps)
end
mh_chain = sample(ordered_multinomial(responses, total_n, 1.5, 4.5), MH(diagm(ones(4) .* [0.1, 0.1, 0.02, 0.02])), 10000)
# NUTS/HMC versions just sit there, even overnight. CPU chart shows it doing SOMETHING, but no idea what
# nuts_chain = sample(ordered_multinomial(responses, total_n, 1.5, 4.5), NUTS(), 1000)
# Check type instability
model = ordered_multinomial(responses, total_n, 1.5, 4.5)
@code_warntype model.f(
model,
Turing.VarInfo(model),
Turing.SamplingContext(
Random.default_rng(), Turing.SampleFromPrior(), Turing.DefaultContext()
),
model.args...,
)
# Multinomial distribution shows up in red, although
# I'm not sure if that's the source of the problem I'm having
```