[Python] Sampling from discrete probability distribution over all natural numbers

I'd like to sample a from the distribution c^a / (-a log (1-c)) where a is a natural number (1 to infinity) and c < 1. How can I do this in Python? I've looked at scipy.stats.rv_discrete, but it seems to require putting in a tuple listing all support values and their corresponding probabilities, and I can't do that for an infinite number of a values.

Continue reading...