Rでべき乗分布(power-law)分布の乱数生成

お仕事で、べき乗分布(power-law)に従った乱数発生をする必要が生じたので、さくっと。

べき乗分布の乱数

検索したら出てきたのがこれ。
Random Number -- from Wolfram MathWorld
べき乗分布に従った乱数 $X$ は、 $y$ が $[0,1]$ の一様乱数、 $n$ がべき乗係数、 $X \in [x_0, x_1]$ のとき、以下となる。
$X = [ (x_1^{n+1} - x_0^{n+1})y + {x_0}^{n+1}]^{\frac{1}{(n+1)}}$

とりあえず、べき乗分布乱数を発生させる関数、rpowerlaw() をほげってみる。

rpowerlaw <- function (n, dist.power = 0, x0 = 0, x1 = 1) {
    sapply(runif(n),
           function (x.unif) (((x1 ^(dist.power + 1) - x0 ^(dist.power + 1))
                               * x.unif + x0 ^(dist.power + 1))
                              ^(1 / (dist.power + 1)))
           )
}

使用例

10000個の、0..1の範囲で分布に従った乱数を生成する場合、とか。

hist(rpowerlaw(10000,dist.power = 7, 0, 1))

f:id:alcuin:20130325042739p:plain

うむ、ろんぐてーるー。

参考

Power-law Distributions

This page is a companion for the SIAM Review paper on power-law distributions in empirical data, written by Aaron Clauset (me), Cosma R. Shalizi and M.E.J. Newman. This page hosts implementations of the methods we describe in the article, including several by authors other than us. Our goal is for the methods to be widely accessible to the community.