Is it possible to model systems with random variables of unknown distribution?

I am new to PyMC and I’ve seen some hierarchical models build from standard distributions that can be found in a statistics book (Bernoulli, Poisson, Exponential, Normal, etc.) Do I get it right that this is always the case in PyMC: one always has to specify a known distribution for all parts of the modelled system? What do people do if some part of the system is poorly approximated by standard distributions?

You can try to define it as a mixture of known distributions or define the dist and logp manually. See for instance the discussion here which could be interesting:

or here for custom dists

If it is difficult to define the likelihood or some of the priors with a closed form formula maybe approximate bayesian computation could be one way:

I do not have experience with ABC but the simulator class seems quite flexible in what it can do. See for instance the code in the first message here (where the op adds random noise to the model that is neither coming from a likelihood or is modelled as a prior):

You might also wanna check things like histogram approximation, bspline_interpolation in pymc_experimental: