A brief answer (sorry I dont have much time to dig deep into it):
-
Yes this is normal sampling from a mixture. There is no good solution in high dimension unfortunately. One possible is to compute the posterior covariance matrix of all mixture component and use that as mass matrix. Because otherwise NUTS will just adapt to the local geometry and explore one local mode. But doing that you might see sampling problem like poor mixing etc.
-
Unfortunately you need to use a for loop as we dont have good support for high dimensional multivariate normal.