You want to get into the habit of using the `pytensor`

functions instead of numpy inside pymc models. Sometimes things work as expected, but it’s not technically supported. In particular, all of those `np.array`

calls with pymc variables will not work as expected.

I guess the only thing missing is `meshgrid`

, which just looks like this:

```
def pt_meshgrid(*xi):
ndim = len(xi)
s0 = (1,) * ndim
outputs = [x.reshape(s0[:i] + (-1,) + s0[i + 1:]) for i, x in enumerate(xi)]
return pt.broadcast_arrays(*outputs)
```

So you should rewrite you code as:

```
with pm.Model() as basic_model:
x = pm.Normal('xval', mu=0, sigma=1)
y = pm.Normal('yval', mu=1, sigma=2)
z = pm.Normal('zval', mu=2, sigma=3)
angle1 = pm.Normal('angle1', mu=0, sigma=np.pi/4)
angle2 = pm.Normal('angle2', mu=0, sigma=np.pi/4)
angle3 = pm.Normal('angle3', mu=0, sigma=np.pi/4)
#Perform some rotations of the x,y,z here
phi = pt.linspace(0,2 * np.pi, 100)
theta = pt.linspace(0, np.pi, 50)
phi, theta = pt_meshgrid(phi, theta)
xc = x*pt.sin(theta) * pt.cos(phi)
yc = y*pt.sin(theta) * pt.sin(phi)
zc = z*pt.cos(theta)
R1 = pt.stack([[pt.cos(angle1), -pt.sin(angle1), 0], [pt.sin(angle1), pt.cos(angle1), 0], [0,0,1]])
R2 = pt.stack([[1,0,0], [0, pt.cos(angle2), -pt.sin(angle2)], [0, pt.sin(angle2), pt.cos(angle2)]])
R3 = pt.stack([[pt.cos(angle3), -pt.sin(angle3), 0], [pt.sin(angle3), pt.cos(angle3), 0], [0,0,1]])
rotated_coords = pt.dot(R3, pt.dot(R2, pt.dot(R1, pt.stack([xc.flatten(), yc.flatten(), zc.flatten()]))))
rotated_coords = rotated_coords.reshape((3, -1))
mu = pm.Deterministic('mu', pt.stack([rotated_coords[0][0], rotated_coords[1][0], rotated_coords[2][0]]))
cov = pt.eye(3) * 0.1
y_hat = pm.MvNormal('y_hat', mu=mu, cov=cov, observed=X)
idata = pm.sample()
```

That meshgrid function would make a good first PR, you’d just have to write a couple tests for it