Indexes of RVs in modeling?

It seems that in any given PyMC3 model, PyMC3 determines the appropriate random variable indices, and it’s hard or impossible for a user to specify such. For example, in:

responses = np.array([ 
          [0,1,2,2,2],
          [0,1,2,1,1],
          [0,1,2,0,0],
          [0,1,2,0,1],
          [0,1,2,1,0]  
    ])
with pm.Model() as model:
    z_student = pm.Normal("z_student", mu=0, sigma=1, shape=(5,3))
    z_question = pm.Normal("z_question",mu=0, sigma=1, shape=(3,5))
    theta = pm.Deterministic("theta", tt.nnet.softmax(z_student - z_question))
    kij = pm.Categorical("kij", p=theta, observed=responses)

The last line kij = pm.Categorical("kij", p=theta, observed=responses) tells PyMC3 that the observed data can be found in a matrix, responses, but it’s not told at all how to retrieve/compute the correct parameters, theta.

It’s implicit that ith row and jth column of responses have some information and should be used in some way in computation of theta. But I’m finding it impossible to tell PyMC3 "retrieve the row and column indices, i & j, of the kth entry in responses, and retrieve the ith row in z_student and the jth column in in z_question for computation of theta.

Could someone explain how these explicit row/column instructions might be given to PyMC3?