重み付き乱数列の生成 (Java)

個人的なコードを書いている中で、重みを付けた乱数列が必要になりました。
具体的には、数が大きいほど発生頻度が高くなるような乱数の生成です。

これについて、どうするのが一般的なのか知らなかったので、調べてみました。
（結論から言うと、一般的な方法がこれ、というのは分かっていません。）

計算式自体は言語を問いませんが、実装例はJavaを使っています。

実行環境

サンプルコードでストリームやラムダを使っているため、Java8以上が必要です。

Java8
- Windows 7 64bit
- JREu101

重み付け無しで乱数列を生成

intの乱数列を取得するのに、以下のコードを実行しました。
IntSupplierのラムダ式を変更すれば、乱数生成方法を変えることができます。

乱数生成関数を用意する必要が無ければ、java.util.Randomクラスのints()メソッドで乱数のIntStreamが得られます。

乱数列（重み付け無し）をストリームで生成

// import java.util.*;
// import java.util.stream.IntStream;

// final int min = 0; // minは省略
final int max = 9; // 実際はもっと大きな数を使う
final long limit = 10000;

final int maxExclusive = max + 1;

// TODO これを実装
IntSupplier weightedRandomGenerator = () -> (int)(Math.random() * maxExclusive);

int[] values = IntStream.generate(weightedRandomGenerator).limit(limit).toArray();

// ここから下は結果を表示
System.out.println(Arrays.toString(values));
for (int i = 0; i < maxExclusive; i++) {
    final int i0 = i;
    final int count = (int)IntStream.of(values).filter(x -> x == i0).count();
    char[] chars = new char[count / 20];
    Arrays.fill(chars, 'o');
    System.out.printf("%2d: %6d %s%n", i, count, String.valueOf(chars));
}

実行結果（重み付け無し）

[3, 2, 0, 6, 3, 7, 1, 6, 1, 3, 9, 5, 8, 6, 1, 2, 7, 4, 7, 9, ... // 省略
 0:   1038 ooooooooooooooooooooooooooooooooooooooooooooooooooo
 1:   1002 oooooooooooooooooooooooooooooooooooooooooooooooooo
 2:   1026 ooooooooooooooooooooooooooooooooooooooooooooooooooo
 3:   1004 oooooooooooooooooooooooooooooooooooooooooooooooooo
 4:   1001 oooooooooooooooooooooooooooooooooooooooooooooooooo
 5:    954 ooooooooooooooooooooooooooooooooooooooooooooooo
 6:    954 ooooooooooooooooooooooooooooooooooooooooooooooo
 7:   1040 oooooooooooooooooooooooooooooooooooooooooooooooooooo
 8:   1006 oooooooooooooooooooooooooooooooooooooooooooooooooo
 9:    975 oooooooooooooooooooooooooooooooooooooooooooooooo

視覚的に分かるように、擬似的なグラフを出力しています。

ほぼ誤差の範囲内で均等にばらけているように見えます。

重み付けのアルゴリズム？を調べる

この辺の知識に乏しいため、どうするのか取っ掛かりが無く、かなりてきとーにググってみたら、Stack Overflowの下記の質問が見つかりました。

random - Weighted randomness in Java - Stack Overflow
- http://stackoverflow.com/questions/6737283/weighted-randomness-in-java

なるほど、「指数分布」(exponential distribution)を使うのですね。

※なお、指数分布の例が多いだけで、必ずしも指数分布でなくても良いみたいです。

指数分布 - Wikipedia
- https://ja.wikipedia.org/wiki/%E6%8C%87%E6%95%B0%E5%88%86%E5%B8%83

さらに検索すると、このサイトを見つけました。
Javaで書いてあるし説明も細かい。ちょうど良さそう。

４．　どんな分布の乱数でも作り出せる－－乱数の変換方法－－（石川宏ナチュラル研究所）
- http://www.ishikawa-lab.com/montecarlo/4shou.html

あとはたくさんありますが、例えばこの辺。

エクセルを用いた指数分布に従う乱数(指数乱数)の発生とヒストグラム
- http://homepage1.nifty.com/gfk/Exp-Rand.htm

モデル化とシミュレーション10
- http://www.junko-k.com/jyoho/simulation/sisuuransuu.htm

これらの情報を元にすると、指数分布に従う指数乱数を生成する数式は

$τ=-\frac{1}{λ}log(1-F(x))$ （ $F(x)$ は一様乱数を生成する関数）

で良さそうです。

あとはこれをコードに変換するだけ...のはず。

実装

数学的な説明は苦手なので省きますが、関数の曲線（グラフ）を書くイメージですね。

二次関数

指数分布では無いですが、ごく単純に $y = x^2$ （ $x >= 0$ ）の二次曲線を書くように関数を定義すれば、それっぽくなります。
これを何と呼ぶのかが分からないので、そのまま期待値が二次曲線になる分布と呼びます。

二次関数による実装

IntSupplier weightedRandomGenerator = () -> {
    double rn = Math.random(); // x
    double biased = rn * rn; // y = x * x
    return (int)(biased * maxExclusive);
};

実行結果（二次曲線による実装）

[0, 5, 0, 1, 6, 0, 9, 6, 1, 0, 0, 0, 0, 1, 6, 8, 9, 8, 3, 7, 5, ... // 省略
 0:   3197 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
 1:   1265 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
 2:   1014 oooooooooooooooooooooooooooooooooooooooooooooooooo
 3:    838 ooooooooooooooooooooooooooooooooooooooooo
 4:    762 oooooooooooooooooooooooooooooooooooooo
 5:    696 oooooooooooooooooooooooooooooooooo
 6:    604 oooooooooooooooooooooooooooooo
 7:    574 oooooooooooooooooooooooooooo
 8:    537 oooooooooooooooooooooooooo
 9:    513 ooooooooooooooooooooooooo

大小が逆ですが、重み付けはそれなりに出来ているっぽいです。

指数乱数 (1)

今度はちゃんと、指数乱数を使って実装します。

指数乱数による実装

final double lambda = 1d;
IntSupplier weightedRandomGenerator = () -> {
    double t = -1d / lambda * Math.log(1d - Math.random());
    return max - (int)t;
};

一部の参考資料の説明では、一様乱数を表す1-F(x)はF(x)に簡略化していますが、Math.random()は0 ≦ x ＜ 1を返し、また、Math.log(0d)は-Infinityになりますので、ゼロにならないように1-F(x)のままにしています。

重み付けは、大きい数値の方が多く出やすくしたいので、max - tとしています。

実行結果（指数乱数による実装）

[7, 8, 9, 8, 9, 8, 9, 7, 6, 9, 8, 6, 8, 8, 9, 8, 8, 9, 8, 8, 6, ... // 省略
 0:      0 
 1:      3 
 2:      8 
 3:     16 
 4:     47 oo
 5:    122 oooooo
 6:    315 ooooooooooooooo
 7:    828 ooooooooooooooooooooooooooooooooooooooooo
 8:   2350 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
 9:   6311 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

ちょっと傾きが極端すぎますが、あとはλの値を調整すれば、求めているものは達成できそうです。

指数乱数 (2)

実は、前項の実装では、少し問題があります。
このコードだと、ゼロ未満の値が稀に（1/1000くらい？）生成されてしまうのです。
理由は単純で、Math.log(x)のxがゼロに近すぎると、値が大きくなりすぎてしまうからです。

そして、この計算では範囲を指定することができません。

私が求めているものはそれほど厳密さを要求しておらず*1、範囲外の値は捨ててしまえば良いので、前項の結果で事足りています。
...が、ここまでやったので、もう少しだけ詰めておきましょう。

と言っても、ここは自力では解決できそうにないので、teratailで質問してみました。

アルゴリズム - ランダム値の重み付けに指数分布を使いつつ範囲を指定したい(45684)｜teratail
- https://teratail.com/questions/45684

質問中では使用している式が違いますが、考え方は理解できました。

一様乱数の範囲を狭めることで、指数乱数自体が0 ≦ x ＜ 1を返すようにすれば、Math.random()のような乱数生成関数と同じように使うことができます。
これで、もう少しだけ精度を高めることができそうです。

ただし、この変形を行うと傾きが小さすぎるので、lambdaを10.0に変えます。

-1d / lambda * Math.log(x)の結果が0 ≦ x ＜ 1になるようなxの範囲は、lambdaが10.0の場合、小数点第11位までの精度にすると0.00004539993～0.99999999999 （0.9999999972741～0.0000000000010）にすると良さそうです。両端の数はやや確率が下がることになりますが、誤差の範囲でしょう。

指数乱数による実装（改良版）

final double lambda = 10d;
final double rmin = 0.00004539993d; // 範囲の下限
final double rmax = 0.99999999999d; // 範囲の上限
final double rw = rmax - rmin; // 範囲の長さ
IntSupplier weightedRandomGenerator = () -> {
    double x = Math.random() * rw + rmin;
    double t = -(1.0d / lambda) * Math.log(x);
    assert 0d <= t && t < 1;
    final int rn = (int)(t * maxExclusive);
    return max - rn;
};

実行結果（指数乱数による実装（改良版））

[9, 8, 9, 9, 9, 8, 8, 9, 8, 9, 9, 9, 7, 9, 9, 8, 9, 9, 8, 8, 9, ... // 省略
 0:      1 
 1:      1 
 2:      5 
 3:     18 
 4:     39 o
 5:    123 oooooo
 6:    325 oooooooooooooooo
 7:    887 oooooooooooooooooooooooooooooooooooooooooooo
 8:   2279 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
 9:   6322 oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

だいたい最初のと同じ傾きになりました。
tも0 ≦ x ＜ 1の範囲に収まることを確認しています。

以上。けっこう時間がかかりました。

数学できないとこういう風に手間取ってしまうので、
基礎くらいはしっかり身に着けておいたほうが良いですね。

（おわり）

*1:少しくらいゆがんだサイコロでもOKな用途。

argius note

プログラミング関連