argius note

プログラミング関連

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

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

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


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


実行環境

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





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

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の下記の質問が見つかりました。


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

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



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

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


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

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

で良さそうです。

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




実装

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


二次関数

指数分布では無いですが、ごく単純にy = x^2x >= 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で質問してみました。

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

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



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

-1d / lambda * Math.log(x)の結果が0 ≦ x < 1になるようなxの範囲は、lambda10.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


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



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

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


(おわり)

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