取りうる値の個数が有限個の任意の離散分布に従う乱数を発生させる「別名法 (alias method)」を Perl で実装してみました。ロジックは下記参考文献に載っていたのそのままで、ソース中のコメントは引用となっています。

■東京大学教養学部統計学教室 (編集), "自然科学の統計学", 東京大学出版会, 1992.


別名法は、例えば「大吉15%、中吉30%、吉30%、凶20%、大凶5%」の割合でランダムにおみくじを出すプログラムを書くときの効率的なアルゴリズムです。

深く考えない実装だと、

$r = rand(1)
if ($r < 0.15) {
  return "大吉";
} elsif ($r < 0.45) {
  return "中吉";
} elsif ($r < 0.75) {
  return "吉";
} elsif ($r < 0.95) {
  return "凶";
} else {
  return "大凶";
}

みたいになるかと思います。このおみくじは5種類だからそれほどではないけど、これが100種類とか1000000種類とかだとどうでしょう。種類数に比例した数の比較演算が必要で時間がかかります。O(N)ですね。（ちょっと工夫して二分探索を使うと O(logN) になります。）

ところが、別名法では準備にほんの少し時間がかかりますが、それ以降は種類数によらずたった1回の比較演算で済みます。O(1)です。すごいですねえ。

どういう仕組みなのかというと……。

って、余裕があるときに「図を用いた分かりやすい説明」を書こうと思います。のんびりお待ちください。

以下、コードと実行例です。

■動作確認スクリプト (aliasmethod.pl)：

#!/usr/bin/perl
use strict;
use warnings;

###### 別名法 (alias method)
### 取りうる値の個数が有限個の任意の離散分布
###   P(X=x_{k}) = p_{k} (1 <= k <= K)
### に従う乱数を発生させるためにウォーカー(Walker, A.J.)が考案した
### 巧妙な方法である。逆関数法と違って、K がいくら大きくても
### 1個の乱数の発生に必要な比較演算が1回だけであるという特徴を持つ。
### ただし、そのために最初に次のように表を用意しておく必要がある。
### (参考文献："自然科学の統計学", 東京大学出版会, 1992)

my @ns = @ARGV;
my $K = @ns; # K
my $sum = 0;
foreach my $i (@ns) {
    $sum += $i;
}
my @p = map {$_ / $sum} @ns; # p_{k}
print "p_{k} :@p\n\n";

### [初期設定] 次の手順に従って a_{k}, v_{k} (1 <= k <= K) を求める

### 1 : v_{k} ← K*p_{k} (1 <= k <= K)
my @v = map {$K * $_} @p; # v_{k}

### 2 : v_{k} >= 1 を満たす k の集合を G,
###     v_{k} < 1 を満たす k の集合を S　とする。
my @G;
my @S;
foreach my $k (0..$K-1) {
    push @G, $k if $v[$k] >= 1;
    push @S, $k if $v[$k] < 1;
}

### 3 : S が空集合でない限り、3.1 から 3.4 までを繰り返し実行する。
my @a; # a_{k}
while (@S and @G) {
    ### 3.1 : S, G から1つずつ任意の要素を選ぶ。それを k, l とする。
    my $k = $S[0];
    my $l = $G[0];

    print "v_{k}: @v\n";
    print "S: @S\n";
    print "G: @G\n";
    print "k=$k, l=$l\n";
    print "\n";

    ### 3.2 : a_{k} ← x_{l}, v_{l} ← v_{l} - (1 - v_{k})
    $a[$k] = $l;
    $v[$l] = $v[$l] - (1 - $v[$k]);

    ### 3.3 : S から k を除く。
    shift @S;

    ### 3.4 : v_{l} < 1 なら l を G から S に移す。
    if ($v[$l] < 1) {
        push @S, (shift @G);
    }
}

print "k p_{k} Kp_{k} v_{k} a_{k}\n";
for (my $i = 0; $i < $K; $i++) {
    print "$i $p[$i] ".($K*$p[$i])." $v[$i] ".(defined $a[$i]?$a[$i]:"*")."\n";
}
print "\n";

### [乱数発生]
sub gen {
    my ($a_r, $v_r) = @_;
    ### 1 : (0, K) 上の一様乱数 V = KU を発生する。
    my $V = rand(@$v_r);
    ### 2 : k ← L V 」+ 1, u ← k - V
    my $k = int($V);
    my $u = $V - $k;
    ### 3 : u <= v_{k} ならば X ← x_{k}, そうでなければ X ← a_{k}
    if ($u <= $v_r->[$k]) {
        return $k;
    } else {
        return $a_r->[$k];
    }
}

### 実験
my $N = 1000000;
my @count;
for (my $i = 0; $i < $N; $i++) {
    my $X = gen(\@a, \@v);
    $count[$X]++;
}
for (my $i = 0; $i < $K; $i++) {
    print "$i ".("*"x($count[$i]/10000))." $count[$i]\n";
}

■実行例：

% aliasmethod.pl 0.12 0.24 0.38 0.16 0.10
p_{k} :0.12 0.24 0.38 0.16 0.1

v_{k}: 0.6 1.2 1.9 0.8 0.5
S: 0 3 4
G: 1 2
k=0, l=1

v_{k}: 0.6 0.8 1.9 0.8 0.5
S: 3 4 1
G: 2
k=3, l=2

v_{k}: 0.6 0.8 1.7 0.8 0.5
S: 4 1
G: 2
k=4, l=2

v_{k}: 0.6 0.8 1.2 0.8 0.5
S: 1
G: 2
k=1, l=2

k p_{k} Kp_{k} v_{k} a_{k}
0 0.12 0.6 0.6 1
1 0.24 1.2 0.8 2
2 0.38 1.9 1 *
3 0.16 0.8 0.8 2
4 0.1 0.5 0.5 2

0 *********** 119632
1 ************************ 240086
2 ************************************* 379950
3 *************** 159991
4 ********** 100341

■短いソース：お作法的にはアレ（特にmap）だけど短い版。

#!/usr/bin/perl
use strict;
use warnings;

my @ns = @ARGV;
my ($a_r, $v_r) = pre(@ns+0, \@ns);

my $N = 1000000;
my @c = do {my @t; map {$t[gen($a_r, $v_r)]++} 1..$N; @t};
print join("", map {"$_ ".("*"x($c[$_]/($N/100)))." $c[$_]\n"} 0..$#ns);

sub pre {
    my ($K, $ns_r) = @_;
    my $sum = do {my $s; map {$s += $_} @$ns_r; $s};
    my (@v, @a, @p, @G, @S);
    map {$_ /= $sum; push @p, $_; push @v, $K * $_} @$ns_r;
    map {$v[$_] >= 1 ? push @G, $_ : push @S, $_} 0..$K-1;
    for (my ($k, $l); $k = shift @S, $l = $a[$k] = $G[0]; ) {
        $v[$l] = $v[$l] - (1 - $v[$k]);
        push @S, shift @G if $v[$l] < 1;
    }
    return \@a, \@v;
}

sub gen {
    my ($a_r, $v_r) = @_;
    my $V = rand(@$v_r);
    my $k = int($V);
    return $V - $k <= $v_r->[$k] ? $k : $a_r->[$k];
}

関連

- unnonouno: 高速な復元抽出の直感的な説明

- discrete_distribution(C++11)

- 非一様乱数 - ORWiki

別名法(alias method)。「取りうる値の個数(n), が有限個の離散分布に従う乱数を, nに無関係なO(1)の速度で生成できる, 簡単で効率的な方法である. ただし, 初期設定にO(n)の時間とO(n)のメモリを必要とする.」

- Alias method - Wikipedia, the free encyclopedia

In computing, the alias method is a family of efficient algorithms for sampling from a discrete probability distribution, due to A. J. Walker. The algorithms typically use O(n log n) or O(n) preprocessing time, after which random values can be drawn from the distribution in O(1) time.

追記140605:
- RSS を読み込んでランダムな文章を生成する[2008-09-07-3]
（別名法を知る前の私の非効率的な実装。確率分布によりランダムにざっくり選ぶ。）