とある偽シャッフルアルゴリズムとその分布

次のようなシャッフルアルゴリズムを考える(簡単のためrand()%Nと表記したが、この部分で0以上N-1未満の一様な整数乱数が生成されると仮定して議論する)。出力されるものは 0, ..., 255 を並び換えたもの(置換)である。

std::vector<int> a(N);
for(int i = 0; i < N; ++i) {
  a[i] = i;
}
for(int i = 0; i < N; ++i) {
  std::swap(a[i], a[rand()%N]);
}

このアルゴリズムは均一ではない。a[i]==jとなる確率は、 i < jのときに高くなり、j <= iのときに低くなる。グラフにすると以下のようになる。

f:id:qnighy:20160605175316p:plain

定理

上のシャッフルアルゴリズムで得られる aについて、 a[i]==jとなる確率は、

${ \begin{cases} \frac{1}{N}\left(1-\frac{1}{N}\right)^{N-1-i} + \frac{1}{N}\left(1-\frac{1}{N}\right)^j - \frac{1}{N}\left(1-\frac{1}{N}\right)^{N-1-i+j} & j \leq i \\ \frac{1}{N}\left(1-\frac{1}{N}\right)^{N-1-i} + \frac{1}{N}\left(1-\frac{1}{N}\right)^j & i < j \end{cases} }$

となる。

証明

2個目のループの本体が ${t}$ 回 ${(0 \leq t \leq N)}$ 実行された時点で a[i]==j となる確率を ${p_{t,i,j}}$ とおく。プログラムをジッと睨むと、以下の式を得る:

${ p_{0,i,j} = \begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases} }$

${ p_{t+1,i,j} = \begin{cases} \frac{1}{N} \sum_{k=0}^{N-1} p_{t,k,j} & i = t \\ \left(1-\frac{1}{N}\right)p_{t,i,j} + \frac{1}{N}p_{t,t,j} & i \neq t \end{cases} }$

さらに、不変条件 ${\sum_{k=0}^{N-1} p_{t,k,j} = 1}$ に注目すると、後者の式は

${ p_{t+1,i,j} = \begin{cases} \frac{1}{N} & i = t \\ \left(1-\frac{1}{N}\right)p_{t,i,j} + \frac{1}{N}p_{t,t,j} & i \neq t \end{cases} }$

となる。このとき、以下が成り立つことを示す。 ${t=N}$ の場合が定理の主張に他ならない。

${ p_{t,i,j} = \begin{cases} \frac{1}{N}\left(1-\frac{1}{N}\right)^{t-1-i} + \frac{1}{N}\left(1-\frac{1}{N}\right)^j - \frac{1}{N}\left(1-\frac{1}{N}\right)^{t-1-i+j} & j \leq i < t \\ \frac{1}{N}\left(1-\frac{1}{N}\right)^{t-1-i} + \frac{1}{N}\left(1-\frac{1}{N}\right)^j & i < j < t \\ \frac{1}{N}\left(1-\frac{1}{N}\right)^j & j < t \leq i \\ \frac{1}{N}\left(1-\frac{1}{N}\right)^{t-1-i} & i < t \leq j \\ \left(1-\frac{1}{N}\right)^t & t \leq i = j \\ 0 & t \leq i, t \leq j, i \neq j \end{cases} }$

${t}$ に関する帰納法で示す。まず、 ${p_{0,i,j}}$ が上記を満たすことはすぐにわかる。

各 ${i,j}$ について ${p_{t,i,j}}$ が上記を満たすと仮定する。このとき、各 ${i,j}$ について ${p_{t+1,i,j}}$ も上記を満たすことを示す。そのために以下のように場合分けをする。

${j \leq i < t}$ のとき。
${j \leq i = t}$ のとき。
${i < j < t}$ のとき。
${i < j = t}$ のとき。
${i < t < j}$ のとき。
${i = t < j}$ のとき。
${j < t < i}$ のとき。
${j = t < i}$ のとき。
${t < i = j}$ のとき。
${t < i, t < j, i \neq j}$ のとき。

各場合分けは簡単な式変形で示せる。

正しいアルゴリズム

例えば、以下のようにする。(別途、rand()%(i+1)の部分をより適切な擬似乱数に置き換える必要がある。)

std::vector<int> a(N);
for(int i = 0; i < N; ++i) {
  a[i] = i;
}
for(int i = 0; i < N; ++i) {
  std::swap(a[i], a[rand()%(i+1)]);
}

この方法の場合、2個目のループ本体が ${t}$ 回実行された時点で a[0], ..., a[t-1] が一様にランダムな置換になっていることが保証できる。

このアルゴリズムはKnuth shuffleとかFisher-Yates shuffleと呼ばれている。

相関に関する補足

この記事では偽シャッフルアルゴリズムについて、要素ごとの確率分布に注目して解析したが、これだけでは不十分な場合もある。例えば、以下のアルゴリズムは明らかに正しくないシャッフルアルゴリズムだが、要素ごとの確率分布は一様である。

std::vector<int> a(N);
int s = rand()%N;
for(int i = 0; i < N; ++i) {
  a[i] = (i+s)%N;
}