Exactly. It’s like BUGS that way. Whereas in PyMC, you have things like Normal("y", mu=intercept + slope * x, sigma=sigma, observed=y) where observed=y is part of the model code for observations whereas something like y = Normal(...) would be used for parameters.
I just did a search through your examples and see you have an example
data = pm.MutableData("data", observed_data[0])
...
pm.Normal("y", mu=mu, sigma=1, observed=data)
which gives you some of what I was looking for in terms of allowing you to plug-and-play different data in the same model (I think). But it still builds observed=data in. So what I was wondering if there’s a way you could write a single PyMC model that you could use for data simulation to do something like simulation based calibration.
The other thing BUGS does which is super-convenient is allow partially instantiated variables (R directly supports this, too, which is costly in compute but useful). So you can do posterior prediction in a regression of y on x by having a partially known y and fully known x. In Stan this all has to get declared up front. We precompile the C++ class for the model with all of the data, parameters and generated quantities declared as such, but we don’t construct an instance until we have data. New data means an entirely new model instance because our models are immutable once constructed.