のんびり読書日記

日々の記録をつらつらと

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で学ぶクラスタ解析

Rで学ぶクラスタ解析

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

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