How to use logp within a model

Is there any way to accomplish a model that works like the following?

with pm.Model():
    a = pm.Uniform('a', lower=-np.pi, upper=np.pi)
    x = pm.VonMises('x', mu=a, kappa=33.3, shape=50)
    y = pm.Deterministic('y', logp({'x[id]': features}), observed=z) 

I have 50 von Mises distributions and I want to calculate the probability in each von Mises for each value of the variable “features”. E.g., for the first von Mises distribution, x[0], I want to compute the probability (y-axis value) for each given x-axis value (given by the variable “features”). I tried to do this via logp but I may be totally wrong in using this approach…

variable legend:
° features = 1000 values ranging from -np.pi to np.pi (if converted from radians to degrees, it would range from 0 to 360 degrees)
° id = 1000 values ranging from 0 to 49 (i.e., each “feature” has an associated id)
° z = 1000 observed values ranging from 0 to 1
° x[id] refers to each of the 50 von Mises distributions