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

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: