Zipf分布に従う乱数の生成方法

Zipf分布といえば，べき乗則でおなじみの分布です（http://en.wikipedia.org/wiki/Zipf%27s_law）．一方，Zipf-MandelbrotはZipf分布を一般化したものだそうです（http://en.wikipedia.org/wiki/Zipf%E2%80%93Mandelbrot_law）．

Zipf-MandelbrotのPDFは

f(x) = c / (x + b)^a

となるわけですが，今回はより単純な形式である，

f(x) = 1 / x

というPDFに従う乱数を生成する方法を説明します．

まず，CDFを求めてみます．ただし，PDFの積分は

∫f(x) dx = ln(x) + C (Cは積分定数)

であるけれど，F(0) = 0なのでCDFは

F(x) = ln(x)

となります．

次に，F(x)の逆関数をもとめます．

y = ln(x)
e^y = x

より，

F^(-1) (x) = e^x

となります．

ここで，逆関数法（http://www.okada.jp.org/RWiki/index.php?%CD%F0%BF%F4Tips%C2%E7%C1%B4#j3f49d96）による乱数生成より，F^(-1) (x)に一様分布に従う乱数を渡せば，Zipf分布に従う乱数が得られます．ただし，実際には[0, 1)の範囲になるように正規化するので，e^x / e = e^(x - 1)を利用します．

Pythonの場合

def zipf():
    return math.e ** (random.random() - 1)

最大値を指定する場合

def zipf(max):
    return math.e ** (random.random() * math.log(max + 1.0)) - 1.0

ちなみに，f(x) = 1 / xという分布は，http://en.wikipedia.org/wiki/Zipf%27s_law#Related_lawsにあるように，英語版のWikipediaの単語分布と非常に近い分布となるらしいです．