べき乗則の乱数
Random numbers generation of the Power Law

■モチベーション

　WEBの世界では観測データがべき乗則に従っている事がしばしばあり、有名なロングテール理論もそれにあたります。例えば、開発しているシステムのテストを行う場合、アクセスパターンがべき乗則に従っている事が分っているならば、べき乗則の乱数でテストデータを作れるわけです。今回はこの様なテストデータを生成するための乱数について考えてみます。

■両対数で直線に

　べき乗則は観測されるデータがパラメータのべき乗に従うことを指します。以下のグラフは y = x ^ -0.7 を両対数でプロットしたもので、この様なデータは両対数上で直線になります。

　log(y) = -0.7 * log(x) ですから -0.7 は傾きを示します。システムのログデータなどから予めこの傾きを把握しておくことが肝要です。

■乱数を作るには

　y = x ^ a から乱数を作るには、x のある範囲における総面積を考え面積 s を一様乱数が与えるパラメータとして x について解けば良さそうです。今回は、x < 1 を取り扱うと s が大きくなるので x ≦ 1 とします。また、積分の都合上 a = -1 とそれ以外で式が異なるので、a = -1 を取り扱いません。使用用途から考えて、a = -0.99 としても実用上問題ないのでその辺はサボる事にします。

■式の展開

　範囲を 1 〜 n として積分し、s をパラメータとする式に変形します。(但し、a ≠ -1)

■R で確認する

　プログラムにするには、s の部分を一様乱数にしてすればOKです。上の式を使って x が 100 万件程度になるような乱数を発生させ、その分布を確認したいと思います。

　a = -0.7 で x = 100万の時の s は、s = (1e6^(-0.7+1) - 1) / (-0.7 + 1) = 206.9858 ですから、乱数は 0 〜 207 の範囲で発生させます。統計ソフトの R に以下のコマンドを入力して100万回試行した分布を確認します。

a <- -0.7
hist( ((a+1)*runif(1e6,0,207)+1)^(1/(a+1)) )

　ヒストグラムに y = x ^ -0.7 を赤色で重ねてみました。大よそ良いみたいです。
Latest update at 2010/4/4