Random Sampling from Custom Distribution

I would like to randomly sample a discrete distribution called the Borel distribution. The Borel distribution is similar to the Poisson distribution, and is defined as:


I have calculated the log likelihood and have implemented it. (I cannot, unfortunately, upload the latex for this as I am a new user.)

Using the class definition of the Poisson distribution as reference, I defined the following class:

import numpy as np
import theano.tensor as tt
from scipy import stats
import warnings

from pymc3.util import get_variable_name
from pymc3.distributions.dist_math import bound, factln, binomln, betaln, logpow, random_choice
from pymc3.distributions.distribution import Discrete, draw_values, generate_samples
from pymc3.distributions.shape_utils import broadcast_distribution_samples
from pymc3.math import tround, sigmoid, logaddexp, logit, log1pexp
from pymc3.theanof import floatX, intX, take_along_axis

class Borel(Discrete):

    def __init__(self, mu, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.mu = mu = tt.as_tensor_variable(floatX(mu))
        self.mode = intX(tt.floor(p))

    def random(self, point=None, size=None):
        #SOME foo
        return None

    def logp(self, value):
        mu = self.mu
        log_prob = bound(
            logpow(mu*value, value-1) - factln(value) - mu*value,
            mu >= 0, value >= 0)
        # Return zero when mu and value are both zero
        return tt.switch(tt.eq(mu, 0) * tt.eq(value, 0),
                         0, log_prob)

How would I calculate the random distribution of the Borel Distribution using PyMC3?

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?