Sample a truncated integer power law in Python?

distribution, numpy, python, random

Solution

AFAIK, neither NumPy nor Scipy defines this distribution for you. However, using SciPy it is easy to define your own discrete distribution function using `scipy.rv_discrete`:

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

def truncated_power_law(a, m):
    x = np.arange(1, m+1, dtype='float')
    pmf = 1/x**a
    pmf /= pmf.sum()
    return stats.rv_discrete(values=(range(1, m+1), pmf))

a, m = 2, 10
d = truncated_power_law(a=a, m=m)

N = 10**4
sample = d.rvs(size=N)

plt.hist(sample, bins=np.arange(m)+0.5)
plt.show()

Problem

What function can I use in Python if I want to sample a truncated integer power law? That is, given two parameters `a` and `m`, generate a random integer `x` in the range `[1,m)` that follows a distribution proportional to `1/x^a`. I've been searching around `numpy.random`, but I haven't found this distribution.

Original source