y_uti のブログ

統計、機械学習、自然言語処理などに興味を持つエンジニアの技術ブログです

Particle filter のリサンプリング方法の比較

Particle filter のリサンプリング方法を比較した下記の論文を読んでみました。そこで、論文中で比較されていた各手法を PHP で実装して試してみたいと思います。

Randal Douc, Olivier Cappé, and Eric Moulines. Comparison of resampling schemes for particle filtering. In 4th International Symposium on Image and Signal Processing and Analysis, 2005.

イメージしやすいように、人工的に作ったでたらめなデータではなく何か実際のデータを使って試したいところです。そう思って探していたところ、Wikipedia に世界各国の人口のリストがあるのを見つけたので、これを使って世界を 100 人の村にしてみることにします。
国の人口順リスト - Wikipedia

まず、ウェブページをダウンロードします。

$ wget http://ja.wikipedia.org/wiki/国の人口順リスト -O world_population.html

「本土と属領合計での順位一覧」のテーブルから国名と人口を抽出します。次のようなコードで比較的簡単に書くことができました。

<?php

function extract_populations($html)
{
    $populations = array();

    $xml = simplexml_load_string($html);
    $targets = $xml->xpath("//table[1]/tr/td[@align='left']");
    foreach ($targets as $t) {
        $c = $t->xpath("a[not(@class='image')][1]");
        if ($c && (string) $c[0] != '世界') {
            $p = $t->xpath("following-sibling::*[1]");
            $populations[(string) $c[0]] = (int) str_replace(',', '', (string) $p[0]);
        }
    }

    return $populations;
}

$populations = extract_populations(file_get_contents('php://stdin'));
foreach ($populations as $c => $p) {
    print "$c: $p\n";
}

実行結果はこんな様子になります。合計は Wikipedia のページの記載と一致しないのですが、気にしないことにします。

$ php extract.php <world_population.html | head -n 20 | column -c 80
中華人民共和国: 1349335152              インドネシア: 239870937
香港: 8053189                           ブラジル: 194946470
マカオ: 557972                          パキスタン: 173593383
インド: 1224514327                      ナイジェリア: 158423182
アメリカ合衆国: 310383948               バングラデシュ: 148692131
プエルトリコ: 3749009                   ロシア: 142958164
グアム: 179896                          日本: 126535920
アメリカ領ヴァージン諸島: 109056        メキシコ: 113423047
アメリカ領サモア: 68420                 フィリピン: 93260798
北マリアナ諸島: 60917                   ベトナム: 87848445
$ php parse.php <world_population.html | awk '{ sum += $NF } END { print sum }'
6881871874

はじめに、最も基本的な multinomial resampling を試してみます。これは、構成比に応じた抽出を必要な回数だけ繰り返す方法です。

<?php

function multinomial_resampling(array $populations, $n)
{
    $samples = array();

    $cumulative = array();
    $total = 0;
    foreach ($populations as $c => $p) {
        $cumulative[$c] = ($total += $p);
    }

    for ($i = 0; $i < $n; ++$i) {
        $r = mt_rand(0, $total - 1);
        foreach ($cumulative as $c => $p) {
            if ($r < $p) {
                @(++$samples[$c]);
                break;
            }
        }
    }

    return $samples;
}

$populations = extract_populations(file_get_contents('php://stdin'));
$samples = multinomial_resampling($populations, 100);
foreach ($samples as $c => $p) {
    print "$c: $p\n";
}

前半のループで $cumulative に $populations の累積値を代入していきます。ループを抜けたときには $total が全体の合計になっています。後半のループでは、0 <= $r < $total の範囲の乱数を一つ生成して国を決めるという処理を $n 回繰り返します。ここは二分探索など高速化の余地がありますが、今回はこだわりません。

このコードを実行して、100 人サンプリングした結果を表示させてみます。なお、この題材では本来は非復元抽出 (同じ人が二回以上選ばれることがない) でなければおかしいのですが、これも主旨とは関係ないので気にしないことにします。

$ php multinomial.php <world_population.html | sort -nrsk2,2 | column -c 80
中華人民共和国: 20      パキスタン: 2           アンゴラ: 1
インド: 18              コンゴ民主共和国: 2     セネガル: 1
ナイジェリア: 4         メキシコ: 2             大韓民国: 1
ロシア: 3               アフガニスタン: 1       マダガスカル: 1
アメリカ合衆国: 3       キルギス: 1             ミャンマー: 1
ブラジル: 2             スリランカ: 1           アルゼンチン: 1
イラン: 2               イギリス: 1             ウガンダ: 1
北朝鮮: 2               ハンガリー: 1           ケニア: 1
フィリピン: 2           ボリビア: 1             カザフスタン: 1
トルコ: 2               ウクライナ: 1           カメルーン: 1
ベトナム: 2             カナダ: 1               マラウイ: 1
インドネシア: 2         ザンビア: 1             ドイツ: 1
イタリア: 2             ハイチ: 1               エジプト: 1
エチオピア: 2           バングラデシュ: 1       ペルー: 1
アルジェリア: 2         サウジアラビア: 1

残念ながら日本人は一人も選ばれませんでした。単純なサンプリングでは、こういうことがよく起きます。全体の 6,881,871,874 人のなかで日本人は 126,535,920 人ですので、比率としては約 1.84% になります。100 人のうち何人が日本人になるか、プログラムの実行を何度も繰り返して様子を見てみると、このようになりました。

$ for i in $(seq 1 1000); do php multinomial.php <world_population.html | grep 日本; done | sort -nk2,2 | uniq -c
    287 日本: 1
    285 日本: 2
    164 日本: 3
     72 日本: 4
     31 日本: 5
      8 日本: 6
      3 日本: 7
      1 日本: 8

合計は 851 なので、残りの 149 回は日本人が選ばれなかったことになります。平均は 1.869 人とほぼ比率どおりの結果になっていますが、やはりかなりのばらつきがあります。用途によっては、このようなばらつきを抑えて、できるだけ元の構成比を維持するようにサンプリングしたいことがあります。

そこで、次に residual resampling を試してみます。これは、先ほどの実行例での日本人のように、本来の構成比に対して少なすぎることがないようにサンプリングする方法です。日本人は全体の 1.84% を占めていたので、100 人を選んだ場合には平均 1.84 人が日本人になります。こういう場合、residual resampling では、小数点以下を切り捨てた 1 人は必ず日本人だと決めてしまいます。各国について同様に最低限の人数を確保してから、残りをランダムに選びます。それではコードを書いてみます。

<?php

function residual_resampling(array $populations, $n)
{
    $samples = array();

    $world_population = array_sum($populations);
    $cumulative = array();
    $total = 0;
    $residual = $n;
    foreach ($populations as $c => $p) {
        $expectation = $p / $world_population * $n;
        if (($floored = floor($expectation)) >= 1) {
            $samples[$c] = $floored;
            $residual -= $floored;
        }
        $cumulative[$c] = ($total += $expectation - $floored);
    }

    for ($i = 0; $i < $residual; ++$i) {
        $r = lcg_value() * $total;
        foreach ($cumulative as $c => $p) {
            if ($r < $p) {
                @(++$samples[$c]);
                break;
            }
        }
    }

    return $samples;
}

前半のループに出現する $expectation が、その国の人が平均何人選ばれるかを表します。この小数点以下を切り捨てた $floored の人数はサンプリングせずに確定させます。$residual で残りの人数を管理し、後半のループでは $residual の分だけサンプリングします。サンプリングの方法は multinomial resampling と同じです。なお lcg_value() は 0 以上 1 未満の乱数を発生させる関数です。線形合同法による乱数なので mt_rand の方が性能がよいらしいですが、浮動小数点数に変換するなど面倒なのでサボってしまいました。

さて、residual resampling の結果を見てみます。先ほどと同様に 1,000 回繰り返してみると、以下のようになりました。

$ for i in $(seq 1 1000); do php residual.php <world_population.html | grep 日本; done | sort -nk2,2 | uniq -c
    433 日本: 1
    368 日本: 2
    140 日本: 3
     49 日本: 4
      6 日本: 5
      3 日本: 6
      1 日本: 7

この方法では日本人が最低 1 人は選ばれるので、合計は 1,000 になります。平均を計算すると 1.840 となり、偶然ながら本来の比率によく一致しました。結果からもわかるように、residual resampling では本来の比率よりも多すぎる選ばれ方になる可能性は依然として残ります。

Stratified resampling*1は、少なすぎることも多すぎることもないようにサンプリングできる方法です。この方法は、まず全体を $n 個の区画に分割し、それぞれの区画から一つずつ要素を選びます。次のようなコードに書くことができます。

<?php

function array_shuffle(array $array)
{
    $result = array();

    $keys = array_keys($array);
    shuffle($keys);
    foreach ($keys as $k) {
        $result[$k] = $array[$k];
    }

    return $result;
}

function stratified_resampling(array $populations, $n)
{
    $samples = array();

    $shuffled = array_shuffle($populations);

    $world_population = array_sum($shuffled);
    $cumulative = array();
    $total = 0;
    foreach ($shuffled as $c => $p) {
        $cumulative[$c] = ($total += $p / $world_population * $n);
    }

    $countries = array_keys($cumulative);
    $c = 0;
    for ($i = 0; $i < $n; ++$i) {
        $r = $i + lcg_value();
        while ($r >= $cumulative[$countries[$c]]) {
            ++$c;
        }
        @(++$samples[$countries[$c]]);
    }

    return $samples;
}

PHP の shuffle 関数はキーを維持してくれないので、キーを維持したままシャッフルする関数を作成しています。この関数の実装にあたっては下記ウェブページを参考にしました。
キーを維持したまま配列をランダムに並べ替える【PHP】 - Programming Magic

Stratified resampling では、最初に対象データをシャッフルしなければサンプリングに偏りが出てしまいます。元々のデータの並び順で末尾の方を見てみると次のようになっています。これらの国々はいずれも人口が小さいので、もし、この並び順で固定してしまうと、すべて最後の一区画にまとめられ、そこから一人しか選ばれないことになります。そうすると、たとえばサンマリノ人とパラオ人が 100 人の村で出会うことは絶対にない、というようなことが起きてしまいます。データの並び順をシャッフルすることで、このような問題を避けることができます。

$ php extract.php <world_population.html | tail -n 5
サンマリノ: 31358
パラオ: 20457
ナウル: 10210
ツバル: 9929
バチカン: 458

それでは結果を見てみます。

$ for i in $(seq 1 1000); do php stratified.php <world_population.html | grep 日本; done | sort -nk2,2 | uniq -c
    254 日本: 1
    655 日本: 2
     89 日本: 3

合計は 998 となり、1,000 回中 2 回は日本人が 1 人も選ばれなかったことになります。最大は 3 人で、4 人以上選ばれることはありません。平均は 1.831 人となりました。日本人は世界人口の 1.84% を占めるので、全体を 100 個の区画に分割すると、シャッフルのされ方によって、日本人はそのうちの 2 区画または 3 区画に含まれます。2 区画の場合は 0 ~ 2 人が選ばれます。3 区画の場合は、真ん中の区画は日本人だけになるので、1 ~ 3 人の日本人が選ばれることになります。

最後に systematic rempling です。これは stratified resampling とほとんど変わりません。Stratified resampling では各区画から 1 人選ぶたびに乱数を生成したのに対して、systematic resampling では最初の区画で一回だけ乱数を生成し、あとは等間隔で選んでいきます。実装は次のようになります。

<?php

function systematic_resampling(array $populations, $n)
{
    $samples = array();

    $shuffled = array_shuffle($populations);

    $world_population = array_sum($shuffled);
    $cumulative = array();
    $total = 0;
    foreach ($shuffled as $c => $p) {
        $cumulative[$c] = ($total += $p / $world_population * $n);
    }

    $countries = array_keys($cumulative);
    $c = 0;
    for ($r = lcg_value(); $r < $n; ++$r) {
        while ($r >= $cumulative[$countries[$c]]) {
            ++$c;
        }
        @(++$samples[$countries[$c]]);
    }

    return $samples;
}

これについても同様に結果を見てみます。

$ for i in $(seq 1 1000); do php systematic.php <world_population.html | grep 日本; done | sort -nk2,2 | uniq -c
    157 日本: 1
    843 日本: 2

Systematic resampling では、最初にオフセットを決めたあとは等間隔に選んでいくので、先ほどの stratified resampling のように、区画の分かれ方によって日本人が 0 人や 3 人になることはありません。平均は 1.843 となりました。

*1:この次に説明する systematic resampling の方法を stratified resampling と呼ぶ場合も多いようです。