Multiple response model, broadcast failure?

Though, I’m running into an error when I try to implement this into my model.

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]  
    ])

students = 5
questions = 5
categories = 3

a = tensor.matrix()
b = tensor.matrix()
elem_sub = a[0,0] - b[0,0], a[0,1] - b[1,0], a[0,2] - b[2,0]  
function = theano.function([a,b], elem_sub)

with pm.Model() as model:
    z_student = pm.Normal("z_student", mu=0, sigma=1, shape=(students,categories))
    z_question = pm.Normal("z_question", mu=0, sigma=1, shape=(categories,questions))
    # Transformed parameter
    theta = pm.Deterministic("theta", tt.nnet.softmax(function(z_student,z_question)))   
    # Likelihood
    kij = pm.Categorical("kij", p=theta, observed=responses)

With the error

TypeError: Bad input argument with name "z_student" to theano function with name "<ipython-input-2-2a16f255dca1>:23" at index 0 (0-based).  
Backtrace when that variable is created:

Expected an array-like object, but found a Variable: maybe you are trying to call a function on a (possibly shared) variable instead of a numeric array?

Needless to say, I’m super confused!