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