I’m not sure what you mean by infinite loop. Can you call pm.draw on your likelihood variable?
Whether to use a CustomDist or not is a modeling question. Right now you are modeling your data as a system of linear equations with correlated errors, centered on fe. Is that your intent? If so, what you have is fine.
You also do not need to use as_tensor_variable in your model. All PyMC objects are already tensor variables.