Random Sampling from Custom Distribution

The random method requires you to implement it yourself as it is general not possible to auto-generate it from the density function. You can call pm.sample(...) on a model with no observation to get prior random sample, but since that essentially based on reject sampling it could be very slow and inefficient (especially in high dimension).

Since it is similar to a Poisson, could you wrap/transform random sample from Poisson in this case?