順番学研究所のご案内

無料ブログはココログ

« 2010年3月 | トップページ | 2010年5月 »

2010年4月

2010年4月 4日 (日)

書籍出版社の数理・番外編―べき分布に従う乱数値のつくりかた

 今回は、ちょっと書籍出版の議論から離れて、テクニカルな話です。数式ばっかりです。

 べき分布に従う乱数を発生させる方法は、前回紹介したようなシミュレーションをやろうと思うと必ず必要になりますが、あまり世間に出回っていないようですね。 kfkfさんからコメント欄で依頼があったように、もしかしたらけっこう、需要があるかもしれませんので、そのやり方について解説します。

結論から

 細かい話が面倒だ、という方のために、まず結論から言いましょう。以下のような式をExcelなり、その他言語にぶち込むと、べき分布に従う乱数が発生します。

 ((M^(-b)-m^(-b))r+m^(-b))^(1/-b)

ただし、

  • M:分布の上限値
  • m:分布の下限値(m>0)
  • b:累積分布関数の指数→パレート指数などと呼ばれているもの(b>0)
  • r:0以上1未満の一様分布乱数(ExcelのRAND関数、BASICならRND関数に対応)

 たとえば、前回のシミュレーションで使った、「上限2,000,000、下限3,000、累積分布関数の指数1.63の分布に従う乱数」なら、以下の式をExcelにぶち込むと、発生できます。

 =((2000000^(-1.63)-3000^(-1.63))*RAND()+3000^(-1.63))^(-1/1.63)

 式の形をみてわかるように、分布の上限値が無限大だと使えません。でも、シミュレーションで使うときは、必ず上限値を設定することになるでしょうから、実用上は問題ないでしょう。

なぜそうなるのか 

 しかし、なぜそうなるのでしょうか。この式は、累積分布関数の性質を利用して、計算しているんです。

 べき分布の確率密度関数p(x)は、一般に以下のように表されます。

 p(x)=Cx^-a  (ただしx>0)

ここで、aは確率密度関数の指数(a=b+1の関係になる)、Cは規格化係数(平たくいえば、分布の下限から上限までの確率の合計が1になるように調整するための係数)です。

 これをxについて積分し、 x<Xの範囲の累積確率を求めるのが、累積分布関数P(X)です。これは次のように表されます。

 P(X)=(C/(-a+1))(X^(-a+1)-m^(-a+1))=(C/(-b))(X^(-b)-m^(-b))

 ここでもmは分布の下限値を、bは累積分布関数の指数を表します。

 このとき、P(X)は分布の下限値からXまでの累積確率なので、0から1の間の値をとります。P(X)は関数なので、0から1の任意のPの値にたいし、ひとつのXが対応することになります。

 このことから、P(X)の逆関数を求めて、0から1の一様分布の乱数に当てはめれば、べき分布に対応する乱数Xが求められることになる、というわけです。こういうのを逆分布関数法、というらしいです。最初に挙げた式は、こうやって求めたものです。

 こんなところで、お役に立ちますでしょうか?

 さて、次回からはまた、書籍出版社の数理に戻ります。

« 2010年3月 | トップページ | 2010年5月 »