2009-09-17

Bayesian Setsを試してみた

Programming

この前YAPC Asia 2009に参加してきたのですが、そこで「はてなブックマークのシステムについて」の発表の中で、「はてブの関連エントリはBayesian Setsを使って計算されている」という話を聞いてBayesian Setsに俄然興味が湧いてきました。Bayesian Setsは以前論文だけ少し読んで、あまりよく分からないまま放置していたのですが、せっかくなのでPerlで作って試してみました。

Bayesian Setsについて詳しくは、以下のリンク先の資料をご参照下さい。

実際に作成したコードは以下の通りです。上記のMatlabのコードを参考にさせていただいています。

#!/usr/bin/perl
#
# Bayesian Sets
#
# Usage:
#  % bayesian_sets.pl input.tsv query1 query2 ..
#
# Reference
# - Paper: http://www.gatsby.ucl.ac.uk/~heller/bsets.pdf
# - Matlab code: http://chasen.org/~daiti-m/dist/bsets/
#

use strict;
use warnings;

use constant {
    MAX_OUTPUT => 20,
};

sub read_vectors {
    my $fh = shift;
    my %vectors;
    while (my $line = <$fh>) {
        chomp $line;
        my @arr = split /\t/, $line;
        my $doc_id = shift @arr;
        my %vector = @arr;
        map { $vector{$_} = 1 } keys %vector;
        $vectors{$doc_id} = \%vector;
    }
    return \%vectors;
}

sub get_parameters {
    my ($vectors, $c) = @_;
    my $average_vector = _average_vector($vectors);
    my (%alpha, %beta);
    while (my ($key, $val) = each %{ $average_vector }) {
        $alpha{$key} = $c * $val;
        $beta{$key}  = $c * (1 - $val);
    }
    return (\%alpha, \%beta);
}

sub calc_similarities {
    my ($vectors, $queries, $alpha, $beta) = @_;
    my %query_vector;
    my $num_query;
    foreach my $query (@{ $queries }) {
        if (!exists $vectors->{$query}) {
            warn "[WARN] '$query' doesn't exist in input documents. Skip it\n";
            next;
        }
        while (my ($key, $val) = each %{ $vectors->{$query} }) {
            $query_vector{$key} += $val;
        }
        $num_query++;
    }
    return if !%query_vector;
    my $length_qvec = scalar(keys %query_vector);
    my %weight_vector;
    foreach my $key (keys %{ $alpha }) {
        my $val = $query_vector{$key} || 0;
        $weight_vector{$key} =
            log(1 + $val / $alpha->{$key})
            - log(1 + ($num_query - $val) / $beta->{$key});
    }
    my %score;
    while (my ($doc_id, $vector) = each %{ $vectors }) {
        $score{$doc_id} = _inner_product(\%weight_vector, $vector);
    }
    return \%score;
}

sub _average_vector {
    my $vectors = shift;
    my %average_vector;
    my $num_vector = 0;
    while (my ($doc_id, $vector) = each %{ $vectors }) {
        while (my ($key, $val) = each %{ $vector }) {
            $average_vector{$key} += $val;
        }
        $num_vector++;
    }
    map { $average_vector{$_} /= $num_vector } keys %average_vector;
    return \%average_vector;
}

sub _inner_product {
    my ($v1, $v2) = @_;
    return 0 if !$v1 || !$v2;
    my @keys = scalar(keys %{ $v1 }) < scalar(keys %{ $v2 }) ?
        keys %{ $v1 } : keys %{ $v2 };
    my $prod = 0;
    foreach my $key (@keys) {
        $prod += $v1->{$key} * $v2->{$key} if $v1->{$key} && $v2->{$key};
    }
    return $prod;
}

sub main {
    my $path = shift @ARGV;
    my @queries = @ARGV;
    if (!$path || !@queries) {
        warn "Usage $0 file query1 query2 ..\n";
        exit 1;
    }
    
    open my $fh, $path or die "cannot open file: $path";
    my $vectors = read_vectors($fh);
    return if !$vectors;
    my ($alpha, $beta) = get_parameters($vectors, 2);  # c=2
    my $score = calc_similarities($vectors, \@queries, $alpha, $beta);
    my $count = 0;
    foreach my $doc_id (sort { $score->{$b} <=> $score->{$a} }
        keys %{ $score }) {
        last if ++$count > MAX_OUTPUT;
        printf "%s\t%.3f\n", $doc_id, $score->{$doc_id};
    }
}

main();

入力データには以下のようなタブ区切りのテキストファイルを用意します。データ例は以前の記事で作成した、Wikipediaの各キーワードをTFIDFで重み付けしたデータです。実際にはvalue部分はいらないのですが、以前の名残りでそのままのフォーマットにしています。

# document_id \t key1 \t val1 \t key2 \t val2 \t ...

龍淵郡  郡      2365    長淵    1857    半島    1555    淵      1460    龍      1313    朝鮮    1137    （/     1126    九      987     串      971     長山 788 リョンヨン      777     モングンポ      777     クミポ  777     朝鮮民主主義人民共和国  777     美浦    733     派      720     教会    714     救面 708 ざんかん        708     チャンサンゴッ  708     金浦    605     史上    595     松川    555     プロテスタント  529     分割    527     地      526 岬 513     里      506     建設    505     苔      500     年      496     南岸    490     当時    487     珪砂    483     黄海    481     再編    468 夢 460     ら      460     北朝鮮  444     宣教    404     1884    397     相      397     ソレ    382     甕      379     白島    360     徐      359 灘 344     教会堂  341     景勝    336     崙      334
神さまこんにちは        ピョンテ        2332    チュンジャ      2332    慶州    1381    映画    1270    人      929     旅      827     チョン・ムソン  777 キム・ボヨン       777     ハン・ミヌ      777     国際    756     ペ・チャンホ    708     脳性    689     ソンギ  639     ミヌ    639     彼      624 祭 615     詩人    614     3       574     妊婦    548     麻痺    539     回      481     念願    468     英      458     39      430     Hello   425 ロード     419     アン    414     家      413     モスクワ        409     電車    399     God     395     16      377     乗車    366     遠足    350 神さま     342     無賃    339     !）     332     さ      331     大人    329     小学生  309     ムービー        298     ベルリン        296     とき 295 実家    289     一      288     原題    286     途方    272     神様    264     韓国    263     公開    259
藤原雅長        年      5208    雅      1285    藤原    1271    近臣    1263    位      1138    従      1118    任      1110    長      957     元年    903 3  861     7       813     院      798     近衛    782     五位上民部権大輔        777     嘉      770     視      764     4       760     下      747 日 747     左      713     隆信    708     月      689     わら    638     叙      629     保      600     8       599     能      595     参議    572 守 556     蔵人    521     久安    515     五      507     少将    506     26      503     顕能    500     在任    460     2       451     ふじ    448 中将       448     19      445     処分    431     久      425     死去    424     任命    422     兼      416     教      406     後期    403     長男 401 この間  377     辞任    377

実際に動かしてみます。入力ドキュメント集合のファイルと、(複数の)クエリを指定します。クエリは入力ドキュメント集合中に含まれている document_id しか使用できません。

% ./bayesian_sets.pl input.tsv トヨタ自動車
トヨタ自動車    140.748
豊田英二        47.205
トヨタ・トヨエース      45.501
5チャンネル     42.303
ポルシェ        40.428
トヨタ・マークX 37.372
日産・マーチ    36.659
フォード・モーター      36.333
三菱自動車工業  33.872
マツダ  33.625
ホンダ・インサイト      32.407
トヨタ・タンドラ        32.004
トヨタ・クルーガー      29.959
ジオ (自動車)   29.293
シボレー・キャバリエ    29.142
トヨタF1        28.830
三菱・ディオン  27.994
いすゞ・エルフ  27.712
ヒュンダイ・ソナタ      27.452
富士重工業      27.111

% ./bayesian_sets.pl input.tsv 原宿警察署
原宿警察署      203.638
滝野川警察署    80.618
牛込警察署      59.605
万世橋警察署    52.241
三田警察署 (東京都)     47.331
城東警察署 (東京都)     43.736
王子警察署      42.404
丸の内警察署    37.594
中野警察署 (東京都)     36.628
広島東警察署    32.846
田無警察署      31.444
豊田幹部交番    31.430
小金井警察署    31.352
尾鷲警察署      30.760
日下部警察署    30.656
山梨県警察      30.290
戸塚町  30.032
阿東幹部交番    29.785
神宮前 (渋谷区) 28.061
倶知安警察署    27.630

ちゃんと関連しているものが取得できていますね。本当は「apple banana」と「apple iPod」のように複数の意味を持つ言葉で、複数のクエリを入れた場合にそのクラスタ内の単語が出てくるのを確かめられるとよかったのですが。

今回使用した入力データのサイズは10万個のドキュメント(wikipediaのキーワード)ですが、実行には25秒ほどかかっています。今回のコードは入力データからのインデックスの作成と、クエリからの類似アイテムの検索の両方を同時に処理させているため動作が遅いですが、インデックスは初めに作成して保存しておけるようにしておけば、そこそこのスピードで検索できるようになると思います。

まだBayesian Setsの肝をあまり理解できていないのですが、

与えられたクエリが含まれていると思われるクラスタをオンデマンドで求め、そのクラスタに含まれている他のアイテムを返す
クエリから求めた重みベクトルと入力ドキュメント集合から構成した行列の積(重みベクトルと全文書の内積)を求めるだけで計算できる

あたりがうれしいところなんでしょうか？一般的な転置インデックスを使用した類似コンテンツを求める方法でも同様のことはできそうな気がするのですが、「オンデマンドにクラスタを求める」というところでゴミが消されるんですかね。ちなみに単純に全文書とのcosine距離で類似したものを取ってきた場合だと、以下のような結果になります。大抵は良好な結果が得られますが、たまにクエリによってはゴミだらけになっちゃってますね。

% ./cos.pl input.tsv トヨタ自動車
トヨタ自動車    1.000
豊田英二        0.625
松井石根        0.610
島崎藤村        0.609
鉄道の歴史 (日本)       0.598
土佐電気鉄道    0.592
日産・マーチ    0.591
新関八洲太郎    0.588
近衛道嗣        0.588
岡村貢  0.587
高田保馬        0.587
8月7日  0.586
仙台東照宮      0.584
3月25日 0.584
飛鳥時代        0.583
三菱自動車工業  0.582
蒲郡競艇場      0.580
富士重工業      0.580
大林組  0.579
ホンダ・インサイト      0.576

% ./cos.pl input.tsv 原宿警察署
原宿警察署      1.000
阿東幹部交番    0.776
周南警察署      0.775
竹原警察署      0.773
高梁警察署      0.771
日下部警察署    0.769
大牟田警察署    0.764
三原警察署      0.758
音戸警察署      0.756
横須賀警察署    0.754
城東警察署 (東京都)     0.753
山口南警察署    0.752
松浦警察署      0.752
王子警察署      0.751
滝野川警察署    0.750
東近江警察署    0.750
つがる警察署    0.748
万世橋警察署    0.743
あわら警察署    0.743
丸亀警察署      0.743

あと、この前のpLSIの記事では結局そのまま放置になってるのですが、せっかくなので今度こそパッケージングしてAlgorithm::BayesianSetsを作ろうと思います。あまり使い道はなさそうですが、ざっくりどんなものか知るのには多少役に立ちそうですし、何より僕自身がCPAN Authorになってみたいので…^^;

追記：weight_vectorを求めるところが一部間違っていたため、コードと実行結果、実行時間を修正しました。ただ実行結果はポイント以外はほとんど違いはないと思います。

2009-08-23

東村アキコ『ママはテンパリスト』1、2巻

本

ママはテンパリスト 1

作者: 東村アキコ
出版社/メーカー: 集英社
発売日: 2008/10
メディア: コミック
購入: 32人クリック: 319回
この商品を含むブログ (255件) を見る

ママはテンパリスト 2

作者: 東村アキコ
出版社/メーカー: 集英社
発売日: 2009/06/19
メディア: コミック
購入: 23人クリック: 70回
この商品を含むブログ (139件) を見る

同じ作者がモーニングで連載している「ひまわりっ〜健一レジェンド」が非常に好きなので、何となく買ってみた。作者の日々の子育ての話が書かれているんだけど、これが非常に面白い！子供ってこんなに面白い反応するものなのか。

あとちょこちょこ知らなかった話も出てきて、それもまた面白い。出産のときに会陰切開をする場合があるらしいけど、それは麻酔なしで切るらしい。痛そうだ…と思ったけど、出産の痛みの方がすごすぎて、切ってもほとんど気にならないぐらいらしい。あと本に出てくる子供は3歳ぐらいになって、普通の食事を食べるようになってからも母乳を飲んでいる。ってそんなにながく母乳って出続けるものなのか。

続きが非常に待ち遠しい漫画だ。

2009-08-16

pLSIを試してみた

Programming

これまでにK-means++とfuzzy c-meansを使用したクラスタリングを試してきましたが、今回はpLSI(probabilistic latent semantic indexing, 潜在的意味インデキシング)によるクラスタリングを試してみようと思います。

pLSIは確率・統計的な枠組みで次元縮約を行う枠組みで、なかなか精度がよいらしく色々な論文で見かけます。Google NewsのレコメンドでもpLSIを使用しており、MapReduceで処理を並列化させて高速に実行しているそうです(論文読んでないので間違っているかも)。また入力ベクトルをあらかじめ重み付けしておく必要がなく、文書であれば単語の頻度をそのまま入力として使用できるのもうれしいところです。

より詳しくは以下のWikipediaのエントリか、書籍をご参照下さい。(書籍は処理結果の表8.4が並びがグチャグチャになってるので注意！)

Probabilistic latent semantic analysis - Wikipedia, the free encyclopedia

Rで学ぶクラスタ解析

作者: 新納浩幸
出版社/メーカー: オーム社
発売日: 2007/11
メディア: 単行本
購入: 9人クリック: 207回
この商品を含むブログ (29件) を見る

また以下の工藤拓さんによるC++実装も参考にさせていただきました。

http://chasen.org/~taku/software/plsi/

今回作成したコードは以下の通りです。Perlでしかも愚直に作ってあるので、大きなデータはちょっと無理です。

#!/usr/bin/perl
#
# pLSI (Probabilistic Latent Semantic Indexing)
# http://en.wikipedia.org/wiki/PLSI
#

use strict;
use warnings;
use Getopt::Long;
use Data::Dumper;

use constant {
    DEFAULT_NUM_ITERATION => 50,
    DEFAULT_BETA          => 0.75,
};

# Format: id\tval1,val2,val3,...
sub read_vectors {
    my $fh = shift;
    my @vectors;
    while (my $line = <$fh>) {
        chomp $line;
        next if !$line;

        my ($id, $vecstr) = split /\t/, $line;
        if (!$id || !$vecstr) {
            warn "Illegal format: $line\n";
        }
        my @vector = split /,/, $vecstr;
        push @vectors, {
            id     => $id,
            vector => \@vector,
        };
    }
    return \@vectors;
}

sub init_random {
    my ($nrow, $ncol) = @_;
    my @results;
    for (my $i = 0; $i < $nrow; $i++) {
        my @arr;
        my $sum = 0;
        for (my $j = 0; $j < $ncol; $j++) {
            my $r = rand();
            push @arr, $r;
            $sum += $r;
        }
        @arr = map { $_ / $sum } @arr;
        push @results, \@arr;
    }
    return \@results;
}

sub expect {
    my ($pdz, $pwz, $pz, $num_doc, $num_word, $num_cluster, $beta) = @_;
    my @scores;
    for (my $id = 0; $id < $num_doc; $id++) {
        for (my $iw = 0; $iw < $num_word; $iw++) {
            my $denominator;
            for (my $iz = 0; $iz < $num_cluster; $iz++) {
                $denominator += $pz->[$iz]
                    * (($pwz->[$iw][$iz] * $pdz->[$id][$iz]) ** $beta);
            }
            for (my $iz = 0; $iz < $num_cluster; $iz++) {
                my $numerator = $pz->[$iz]
                    * (($pwz->[$iw][$iz] * $pdz->[$id][$iz]) ** $beta);
                $scores[$id][$iw][$iz] = $denominator ?
                    $numerator / $denominator : 0;
            }
        }
    }
    return \@scores;
}

sub maximize {
    my ($scores, $documents, $pdz, $pwz, $pz,
        $num_doc, $num_word, $num_cluster) = @_;
    my (@pdz_new, @pwz_new, @pz_new);

    my @denominators;
    my $sum_denom;
    for (my $iz = 0; $iz < $num_cluster; $iz++) {
        my $denominator;
        for (my $id = 0; $id < $num_doc; $id++) {
            for (my $iw = 0; $iw < $num_word; $iw++) {
                $denominator += $documents->[$id][$iw] * $scores->[$id][$iw][$iz];
            }
        }
        push @denominators, $denominator;

        # update pdz
        for (my $id = 0; $id < $num_doc; $id++) {
            my $numerator;
            for (my $iw = 0; $iw < $num_word; $iw++) {
                $numerator += $documents->[$id][$iw] * $scores->[$id][$iw][$iz];
            }
            $pdz_new[$id][$iz] = $numerator / $denominator;
        }
        # update pwz
        for (my $iw = 0; $iw < $num_word; $iw++) {
            my $numerator;
            for (my $id = 0; $id < $num_doc; $id++) {
                $numerator += $documents->[$id][$iw] * $scores->[$id][$iw][$iz];
            }
            $pwz_new[$iw][$iz] = $numerator / $denominator;
        }

    }
    # update pz
    my $denom_sum;
    map { $denom_sum += $_ } @denominators;
    for (my $iz = 0; $iz < $num_cluster; $iz++) {
        $pz_new[$iz] = $denominators[$iz] / $denom_sum;
    }
    return (\@pdz_new, \@pwz_new, \@pz_new);
}

sub plsi {
    my ($documents, $num_cluster, $num_iter, $beta) = @_;
    my $num_doc = scalar @{ $documents };
    my $num_word = scalar @{ $documents->[0] };
    my $pdz = init_random($num_doc, $num_cluster);
    my $pwz = init_random($num_word, $num_cluster);
    my $pz = [ map { 1 / $num_cluster } (0 .. $num_cluster-1) ];

    for (my $i = 0; $i < $num_iter; $i++) {
        my $scores = expect(
            $pdz, $pwz, $pz, $num_doc, $num_word, $num_cluster, $beta);
#        print "pdz: ", Dumper($pdz);
#        print "pwz: ", Dumper($pwz);
#        print "pz: ", Dumper($pz);
#        print Dumper($scores);
        ($pdz, $pwz, $pz) = maximize($scores, $documents, $pdz, $pwz, $pz,
                                     $num_doc, $num_word, $num_cluster);
    }
    return ($pdz, $pwz, $pz);
}

sub main {
    my %option;
    GetOptions(
        \%option,
        'number=i',
        'beta=f',
        'iter=i',
    );
    my $path = shift @ARGV;
    if (!$path || !$option{number} || $option{number} <= 0) {
        print "Usage: plsi.pl -n ncluster [-b beta | -i niter] inputfile\n";
        exit 1;
    }
    my $num_cluster = $option{number};
    my $num_iter = exists $option{iter} ? $option{iter} : DEFAULT_NUM_ITERATION;
    my $beta     = exists $option{beta} ? $option{beta} : DEFAULT_BETA;
    die 'iter must be >0' if $num_iter <= 0;
    die 'beta must be >0' if $beta <= 0;

    open my $fh, "<$path" or die "cannot open: $path";
    my $vectors = read_vectors($fh);
    my @documents;
    foreach my $v (@{ $vectors }) {
        push @documents, $v->{vector};
    }
    if ($num_cluster <= 0 || $num_cluster > scalar(@documents)) {
        die 'number of centers must be >0 and <=number_of_vectors';
    }

    my ($pdz, $pwz, $pz) = plsi(\@documents, $num_cluster, $num_iter, $beta);

    for (my $i = 0; $i < scalar @{ $pdz }; $i++) {
        print $vectors->[$i]{id}."\t";
        print join ' ', map { sprintf "%.4f", $_ } @{ $pdz->[$i] };
        print "\n";
    }
}

main();

実際に動かしてみます。実行結果はP(d|z)だけ表示しており、これをクラスタリングの結果と捉えることができます。

% cat input.txt
1       2,1,0,3,2
2       3,3,1,1,2
3       0,1,2,1,1
4       2,2,3,1,2

% ./plsi.pl -n 2 input.txt
1       0.0376 0.3015
2       0.1847 0.3372
3       0.3211 0.1026
4       0.4566 0.2587

入力文書1, 2はクラスタ2に主に属し、文書3, 4はクラスタ1に属していることが分かります。書籍の説明ではEMの各ステップでP(d|z),P(w|z)の合計値を1.0に正規化してはいなかったのですが、C++実装では正規化されてました。ここは揃えた方が気持ちいい、程度の問題なんでしょうか？とりあえず今回作ったプログラムでは各ステップでの正規化はしていません。

もう少しわかりやすいデータでも試してみると、

% cat input2.txt 
1       1,1,1,0,0,0
2       1,0,1,0,0,0
3       1,1,0,0,0,0
4       0,0,0,1,1,1
5       0,0,0,1,0,1
6       0,0,0,1,1,0

./plsi.pl -n 2 /input2.txt 
1       0.0000 0.4286
2       0.0000 0.2857
3       0.0000 0.2857
4       0.4286 0.0000
5       0.2857 0.0000
6       0.2857 0.0000

文書1,2,3がクラスタ2、文書4,5,6がクラスタ1にきっちり分かれました。ちゃんと動いてそうですね。

またbetaパラメータはEMのスムージング用のパラメータで、デフォルトは0.75になっていますが、これを低い値にするとよりなめらかな分布になるそうです(TemperedEMというらしい)。詳しくは調べていないので、パラメータの説明はC++実装のREADMEを参照してください。

% ./plsi.pl -n 2 input.txt
1       0.0376 0.3015
2       0.1847 0.3372
3       0.3211 0.1026
4       0.4566 0.2587

./plsi.pl -n 2 -b 1.0 input.txt
1       0.0000 0.5428
2       0.2168 0.4099
3       0.2731 0.0008
4       0.5101 0.0465

./plsi.pl -n 2 -b 0.65 input.txt
1       0.2412 0.2435
2       0.3027 0.3033
3       0.1522 0.1509
4       0.3038 0.3023

詳しく調べないで何回か試しただけなのですが、初期値によるブレが大きくて、そんなにいいものではないのかな感じたのですが、どうなんでしょう？10回ぐらい同じ処理をして、その平均を取るとかにした方がいいのかな…。

次はlda(latent Dirichlet allocation, 潜在的ディリクレ配分法)を試すぞ！といいたいところですが、説明読んでもよく分からない…orz なのでとりあえずpLSIをC++で書き直して、bayonにでも組み込んでみようかと思います。

それか、CPANにpLSIのパッケージがなかったので、ちゃんと作り直してAlgorithm::PLSIとかでパッケージングしたらうれしい人いるかな？結局pure Perlじゃ大きなデータは厳しいので、使い物にならないかもしれないけど^^;

2009-08-02

鈴木ともこ『山登りはじめました』

本

山登りはじめましためざせ!富士山編

作者: 鈴木ともこ
出版社/メーカー: メディアファクトリー
発売日: 2009/06/17
メディア: 単行本（ソフトカバー）
購入: 6人クリック: 123回
この商品を含むブログ (41件) を見る

なんとなくアウトドアには憧れがあるんです。ただ大変そうなイメージがあるのと、本気でやろうと思うとシューズやらバッグやらレインウェアやら…とすごく金がかかりそうな気がして、なかなか手がでません。あと近頃は全く運動してないので、体力も心配。

この本は元々は運動が苦手でインドア派だった作者が、いろんな山を登るようになって最終的に富士山に登るまでの話が漫画で書かれています。写真も掲載されていて、綺麗な景色も楽しめる。行く先々での食事もすごくおいしそうで、自分も行きたくなってきます。自分のペースでじっくりのんびり行くのが楽しむコツなんですかね。

この本にも書いてありましたが、日本人に生まれたからには一度は富士山に登ってみたい！。でも富士山って7,8月しか登れないんですね。今年はちょっと無理そうなので、来年登るのを目標にして、この本に書いてある場所に行ってみようかなぁと妄想中。まずは高尾山か鎌倉アルプスですかね。

岳 6 (ビッグコミックス)

作者: 石塚真一
出版社/メーカー: 小学館
発売日: 2008/01/30
メディア: コミック
購入: 1人クリック: 39回
この商品を含むブログ (71件) を見る

そうそう、山の本といえば「岳」。今回この本を読んでて、岳に出てきたフレーズを思い出した。

遠くで見たときずっとキレイになるんだよ
頂上まで登った山は

岳は本格的な山登りなので、自分も登りたいという気持ちにはなかなかならなかったけど、今回の本はもっと気軽に登れそうな気がしてすごくよかった。

2009-07-24

STLのvectorとpriority_queueのソート用比較関数は不等号が逆

Programming

この前自分のソースを読んでいたら、両方とも降順にソートするために作った比較関数なのに何故か不等号が逆になっていて、「うわ、ひどいバグ作っちゃった？！」って慌ててテストしたら問題なし。調べてみると、STLのvectorとpriority_queueのソート用比較関数は何故か不等号が逆になっちゃうんですね。

#include <cstdlib>
#include <vector>
#include <queue>
#include <iostream>

struct StructGreater {
  bool operator() (int a, int b) {
    return a < b;
  }
};

bool funcGreater(int a, int b) {
  return a > b;
}

int main() {
  std::priority_queue<int, std::vector<int>, StructGreater> q;
  std::vector<int> v;
  size_t max = 10;
  for (size_t i = 0; i < max; i++) {
    int n = rand() % 100;
    q.push(n);
    v.push_back(n);
  }
  sort(v.begin(), v.end(), funcGreater);

  for (size_t i = 0; i < max; i++) {
    std::cout << "v:" << v[i] << "\tq:" << q.top() << std::endl;
    q.pop();
  }
  return 0;
}

ランダムに生成した整数列をvectorとpriority_queueに入れて、ソートして出力してみます。StructGreaterとfuncGreaterで不等号は逆なのですが、結果はちゃんと両方とも降順にソートされています。

% g++ sort.cc
% ./a.out
vector:93       queue:93
vector:92       queue:92
vector:86       queue:86
vector:86       queue:86
vector:83       queue:83
vector:77       queue:77
vector:49       queue:49
vector:35       queue:35
vector:21       queue:21
vector:15       queue:15

知ってる人にとっては当たり前なのかもしれないけど、念のためメモメモ。

2009-07-24

Fuzzy c-meansを試してみた

Programming

K-meansは各入力ドキュメントがただ1つのクラスタにのみ属するハードクラスタリング手法ですが、fuzzy c-meansは所属度を持って複数のクラスタへの所属を許すソフトクラスタリング手法です。K-meansは以前に作りましたので、今回はfuzzy c-meansを試したいと思います。

fuzzy c-meansの詳しいアルゴリズムは下のページを参照してください。

Cluster analysis - Wikipedia, the free encyclopedia

とりあえずPerlで実装しました。ソースの半分くらいは、以前作ったK-means++からのコピペです^^;

#!/usr/bin/perl
#
# Fuzzy c-means clustering
# http://en.wikipedia.org/wiki/Cluster_Analysis#Fuzzy_c-means_clustering
#

use strict;
use warnings;
use Data::Dumper;
use List::Util qw(shuffle);
use List::MoreUtils qw(any);

sub choose_random_centers {
    my ($vectors, $num_center) = @_;
    my @ids = keys %{ $vectors };
    @ids = shuffle @ids;
    my @centers = map { $vectors->{$_} } @ids[0 .. $num_center-1];
    return \@centers;
}

sub choose_smart_centers {
    my ($vectors, $num_center) = @_;
    my $cur_potential = 0;
    my @centers;

    # choose one random center
    my $vector = (shuffle values %{ $vectors })[0];
    push @centers, $vector;
    my %closest_dist;
    foreach my $id (keys %{ $vectors }) {
        $closest_dist{$id} = squared_euclid_distance($vectors->{$id}, $vector);
        $cur_potential += $closest_dist{$id};
    }

    # choose each center
    for (my $i = 1; $i < $num_center; $i++) {
        my $randval = rand() * $cur_potential;
        my $center_id;
        foreach my $id (keys %{ $vectors }) {
            $center_id = $id;
            last if $randval <= $closest_dist{$id};
            $randval -= $closest_dist{$id};
        }
        my $new_potential = 0;
        foreach my $id (keys %{ $vectors }) {
            my $dist = squared_euclid_distance(
                $vectors->{$id}, $vectors->{$center_id});
            $closest_dist{$id} = $dist if $dist < $closest_dist{$id};
            $new_potential += $closest_dist{$id};
        }
        push @centers, $vectors->{$center_id};
        $cur_potential = $new_potential;
    }
    return \@centers;
}

# Format: id\tval1,val2,val3,...
sub read_vectors {
    my $fh = shift;
    my %vectors;
    while (my $line = <$fh>) {
        chomp $line;
        next if !$line;

        my ($id, $vecstr) = split /\t/, $line;
        if (!$id || !$vecstr) {
            warn "Illegal format: $line\n";
        }
        my @vector = split /,/, $vecstr;
        $vectors{$id} = \@vector;
    }
    return \%vectors;
}

sub squared_euclid_distance {
    my ($vec1, $vec2) = @_;
    my $size = scalar @{ $vec1 };
    my $dist = 0;
    for (my $i = 0; $i < $size; $i++) {
        $dist += ($vec1->[$i] - $vec2->[$i]) * ($vec1->[$i] - $vec2->[$i]);
    }
    return $dist;
}

sub get_distances {
    my ($vectors, $centers) = @_;
    my %distances;
    foreach my $id (keys %{ $vectors }) {
        foreach my $center (@{ $centers }) {
            my $dist = squared_euclid_distance($vectors->{$id}, $center);
            push @{ $distances{$id} }, $dist;
        }
    }
    return \%distances;
}

sub calc_weights {
    my ($vectors, $centers, $m) = @_;
    my $distances = get_distances($vectors, $centers);
    return if !$distances;

    my $num_center = scalar @{ $centers };
    my %weights;
    foreach my $id (keys %{ $distances }) {
        if (any { $_ == 0 } @{ $distances->{$id} }) {
            foreach my $dist (@{ $distances->{$id} }) {
                push @{ $weights{$id} }, $dist == 0 ? 1 : 0;
            }
        }
        else {
            for (my $i = 0; $i < $num_center; $i++) {
                my $weight;
                for (my $j = 0; $j < $num_center; $j++) {
                    my $x = $distances->{$id}[$i] / $distances->{$id}[$j];
                    $weight += $x * $x;
                }
                $weight **= (-1) / ($m - 1);
                push @{ $weights{$id} }, $weight;
            }
        }
    }
    return \%weights;
}

sub calc_centers {
    my ($vectors, $weights, $num_center, $m) = @_;
    my @centers;

    # initialize centers
    my $veclen = scalar @{ (values %{ $vectors })[0] };
    for (my $i = 0; $i < $num_center; $i++) {
        my @arr;
        for (my $j = 0; $j < $veclen; $j++) {
            push @arr, 0;
        }
        push @centers, \@arr;
    }

    # sum of weights
    my @weight_sums;
    for (my $i = 0; $i < $num_center; $i++) {
        push @weight_sums, 0;
    }
    foreach my $id (keys %{ $weights} ) {
        for (my $i = 0; $i < $num_center; $i++) {
            $weight_sums[$i] += $weights->{$id}[$i] ** 2;
        }
    }
    
    # calc center position
    foreach my $id (keys %{ $vectors }) {
        for (my $i = 0; $i < $num_center; $i++) {
            for (my $j = 0; $j < $veclen; $j++) {
                $centers[$i][$j] += $weights->{$id}[$i] ** 2
                    * $vectors->{$id}[$j] / $weight_sums[$i];
            }
        }
    }
    return \@centers;
}

sub fuzzy_cmeans {
    my ($vectors, $num_center, $num_iter, $m) = @_;
 
    # choose initial centers;
    #my $centers = choose_random_centers($vectors, $num_center);
    my $centers = choose_smart_centers($vectors, $num_center);

    my $weights;
    for (my $i = 0; $i < $num_iter; $i++) {
        #warn "Loop: $i\n";
        $weights = calc_weights($vectors, $centers, $m);
        $centers = calc_centers($vectors, $weights, $num_center, $m);
    }
    return $weights;
}

sub main {
    my ($path, $num_center) = @ARGV;
    if (!$path || !$num_center) {
        print "Usage: fuzzy_cmeans.pl inputfile num_center\n";
        exit 1;
    }
    open my $fh, "<$path" or die "cannot open: $path";
    my $vectors = read_vectors($fh);
    my $num_vectors = scalar keys %{ $vectors };
    if ($num_center <= 0 || $num_center > $num_vectors) {
        die 'number of centers must be >0 and <=number_of_vectors';
    }
    
    my $num_iter = 10;
    my $m = 2;
    my $weights = fuzzy_cmeans($vectors, $num_center, $num_iter, $m);
    print Dumper($weights);
}

main()

__END__

通常のfuzzy c-meansは初めの中心の選択はランダムに行うのですが、上記のソースではK-means++と同じ手法を使用して、なるべく似ていないものを選択するようにしています。

d1      1,2
d2      2,2
d3      4,2
d4      1,1
d5      2,1
d6      4,1

上記のデータを使用して、実際にクラスタリングを実行してみます。

% ./fuzzy_cmeans.pl input.txt 2
$VAR1 = {
          'd3' => [
                    '0.00146446087457261',
                    '0.998535539125427'
                  ],
          'd1' => [
                    '0.99714951936669',
                    '0.00285048063330992'
                  ],
          'd2' => [
                    '0.986032261281484',
                    '0.013967738718516'
                  ],
          'd5' => [
                    '0.986032261281475',
                    '0.0139677387185249'
                  ],
          'd4' => [
                    '0.997149519366688',
                    '0.00285048063331179'
                  ],
          'd6' => [
                    '0.00146446087457242',
                    '0.998535539125428'
                  ]
        };

整形して出力してないので見づらいのですが、ポイントが高いほどそのクラスタに適合していると考えられるので、d1,d2,d4,d5がクラスタ1、d3,d6がクラスタ2と大体分けることができます。

また各ドキュメントに対し、最も適合したクラスタ以外のクラスタへの所属度も出力できています。今回のデータではどちらか一方のクラスタへの所属度が大幅に大きくなってしまいましたが、データによっては同じくらいずつの所属度を複数のクラスタに持つこともあります。データは必ずしも1つの属性のみを持つのではなく、複数の属性を持つこともあるので、fuzzy c-meansのように複数のクラスタへの所属度を表示できるのはやはり便利ですね。

Rで学ぶクラスタ解析

作者: 新納浩幸
出版社/メーカー: オーム社
発売日: 2007/11
メディア: 単行本
購入: 9人クリック: 207回
この商品を含むブログ (29件) を見る

ちなみに僕は「Rで学ぶクラスタ解析」でfuzzy c-meansを調べました。この本は説明も読みやすく、数値での実際の計算例も掲載されているため、非常にわかりやすいです。

ただ、今回のfuzzy c-meansは数式で一部誤植があるような？他にもEMアルゴリズムの数式も多分間違っているところがあって、サポートページとか誤植の一覧ページがあると助かるのですが…。

のんびり読書日記

日々の記録をつらつらと

Bayesian Setsを試してみた

東村アキコ『ママはテンパリスト』1、2巻

最近読んだ本

pLSIを試してみた

鈴木ともこ『山登りはじめました』

STLのvectorとpriority_queueのソート用比較関数は不等号が逆

Fuzzy c-meansを試してみた