前回のプログラム(の断片)では、データを乱数で適当にクラスタ分けしました。これは結局、データの分布範囲内で適当に初期の中心点を選んだ事に他なりません。初期値依存、とは、この最初の中心点の選び方でクラスタリング結果が大きく変わってしまう、という事です。
例えば、下の図の2つの円(グレーでハッチされた円)内に、同じデータ密度で同数のデータが一様に分布しているとします。
これを2つのクラスタに分けろ！と言われれば、普通右と左の円で分けるでしょうね。
で、ランダムに選んだ初期中心点が、たまたま図の赤・青丸のところに来たとしたら、それぞれのデータはいちばん近い中心点のクラスタに分けられますから下の図のように直観と同じ結果が得られるでしょう。
でもこれが下の図のように、2つの円の中心を通る線(破線)から同じ距離で、上下に分布してしまったら
こんな結果になる事でしょう。
データ間の距離でクラスタ分けしたらどちらの結果も「正しい」ものなのでしょうが、初期中心点の選び方で結果は大きく異なります。
それに後者の分け方は、なんでわざわざそんな分け方を・・・な感じで、なんだか不自然です。
こういう「不自然な」分け方にならないよう、初期値を上手に選ぶ方法があるそうです。
それがこちらのページで紹介されているkmeans++です。
元々の提案者はこちらの方だそうで論文も公開されています。
まぁ、論文読んでもダビット君が何言ってるのかは分からないがっ！
何をしろと言ってるのかは分かった気がする(?)
論文でやれっと言ってることは要するに・・・
初期の中心点を選ぶとき、各データに最も近い中心点までの距離を
として、以下の確率
を最大にする点を選びなさい
そう言ってる気がする。
確率とか出てくるとなんだか分からなくなるのですが（涙）
でも上の「確率」って、分母は全てのデータについてDの和を取ってるだけなので、どの点でも値は変らないですよね？という事は分子が最大になる点を選べば良い訳？
なのでさらに要するに
「最も近い中心点までの距離」が最大になるデータ点を
次の中心点として選びなさい！
そう言われてる気がする・・・
やり方としては、初期中心点をn個決めなきゃいけない (データをn個のクラスタに分ける) 場合
1個目の中心点は、データの中からまぁテキトーに決める。
以降、それまで決まっている中心点のうち、一番近い中心点までの距離が最大になるデータ点を次の中心点に選ぶ。
こうして全部でn個中心点が決まるまで続ける。
但し中心点は重複させない事！
4個目と6個目の中心点が同じ点だったよ。あはは♪なんてこと無いように。
例えば下のA,B,C,D 4個の点から中心点を順に選ぶとします。
点が4個しかないのにクラスタリングとかいみねくね？とか言わない。
1個目の中心点は適当にDと決めます。
続く2個目の中心点ですが、中心点は重複させない、あはは♪なので残るA,B,Cの中から選びます。A,B,Cから一番近い中心点、は今一個しか中心点がないのでDです。なので最近中心点であるDまでの距離が一番大きいのはA。なので2個目はA
3個目は残るB、Cから選びます。Bに一番近い中心点はA、CはD。距離|BA|と|CD|で最大なのは|BA|。なので3個目はB。
こんな感じで中心点を選び、それを初期状態にして後はクラスタリング＆中心点更新、を普通のkmenas法と同じやり方で繰り返します。
ここの部分のコードを書くとしたら下みたいな感じでしょうか？
前回の構造体 clustered_data *cdata の初期化の(前回は乱数でクラスタ分けしてた)ところで・・・
gsl_vector *c = gsl_vector_alloc (cdata->data->size2);
// 各データが既に中心点に使われたかを表すフラグ
// 中心点の重複は許されない。あはは♪
int *data_is_used
= (int *) malloc (cdata->data->size1 * sizeof (int));
for (i = 0; i < cdata->data->size1; i++) data_is_used[i] = FALSE;
// 最初の中心点を選ぶ
gsl_rng *r = gsl_rng_alloc (gsl_rng_default);
int firstcentroid = gsl_rng_uniform_int (r, cdata->data->size1);
// 乱数で選ばれたfirstcentroid番目のデータ点を最初の中心点に
gsl_matrix_get_row (c, cdata->data, firstcentroid);
gsl_matrix_set_row (cdata->centroids, 0, c);
data_is_used[firstcentroid] = TRUE;
// 残りの初期中心点 cdata->nclusters - 1個を求める
for (k = 1; k < cdata->nclusters; k++) {
int nextcentroid;
double rmax = 0.;
// 各データに最も近い中心点までの距離を求める。
for (i = 0; i < cdata->data->size1; i++) {
double r;
// i番目データに最も近い中心点までの距離を格納
double rmin = GSL_DBL_MAX;
// i番目のデータが既に中心点に使われていたら除外
if (data_is_used[i]) continue;
// i番目のデータを取り出し
gsl_vector *xi = gsl_vector_alloc (cdata->data->size2);
gsl_matrix_get_row (xi, cdata->data, i);
// これまで選ばれた中心点について
for (j = 0; j < k; j++) {
// j番目の中心点を取り出し
gsl_matrix_get_row (c, cdata->centroids, j);
// i番目データとj番目中心点の距離を計算
r = eucledian_distance (xi, c, NULL);
// 今まででi番目データに最も近い中心点までの距離
if (r < rmin) rmin = r;
}
// 上で求まった、i番目データに最も近い中心点までの距離について
if (rmax < rmin) {
// 今までのデータの最近中心点距離の最大値
rmax = rmin;
// rmaxを叩き出したデータ点番号
nextcentroid = i;
}
gsl_vector_free (xi);
}
// 最近中心点距離が最大になるnextcentroid番目のデータを
// 次の初期中心点に設定
gsl_vector *nextc = gsl_vector_alloc (cdata->data->size2);
gsl_matrix_get_row (nextc, cdata->data, nextcentroid);
gsl_matrix_set_row (cdata->centroids, k, nextc);
gsl_vector_free (nextc);
// 使用済み
data_is_used[nextcentroid] = TRUE;
}