Thanks for such an encouraging and in-depth answer!
I fully agree with your way of thinking on the indexer - I get how we’re essentially using numpy advanced indexing to associate the model parameters with to the appropriate observed data.
Your first answer matches the approach used in R in the solutions (page 4, first code block).
Your answer also helps me with my understanding on how to think of the shapes and dimensions in the model (I asked a question about this a few weeks ago).