id:audioswitchさんのX-meansをPHPで実装してみました
id:audioswitchさんのX-meansをPHPで実装してみました。
普段、処理は全部PHPで書いているのでPHPで書き起こしました。
今回はコードも公開して頂いていたので、使わせてもらうだけではなくもしかすると需要があるかもしれないので公開させて頂きました次第であります。
コピーするだけでも能がないので、何かないものかと。おまけでGDでプロットして確認する関数も作ったので一緒に。
あと、サンプルで12ポイントだけプロットしてみるコードもご一緒にどうぞ。
<?php //------------------------------------- // x-means法でクラスタリング public function xmeans_clustering ( array $data, $dimension_number, $sample_number, $init_cluster_number) { $label =array(); $next_label =0; $old_cluster_number =$init_cluster_number; $new_cluster_number =0; // 最初のkmeansを行う． list($old_centroid, $old_label) =kmeans_clustering( $data, $dimension_number, $sample_number, $old_cluster_number); // クラスタ数に変化がなくなるまで繰り返す． do { // 新しいクラスタを初期化する． $new_label =array_diff($old_label,array()); $new_centroid =array_diff($old_centroid,array()); $old_cluster_number; $next_label = 0; // 各クラスタを２つに分割し，評価していく． for ($n=0; $n<$old_cluster_number; $n++) { // n番目のクラスタのデータ数を数える． $parent_sample_number = 0; for ($i=0; $i<$sample_number; $i++) { if ($old_label[$i] == $n) { $parent_sample_number++; } } if ($parent_sample_number <= 1) { continue; } // 分散を求める． $parent_variance =0; for ($i=0; $i<$sample_number; $i++) { $distance =0; if ($old_label[$i] == $n) { for ($j=0; $j<$dimension_number; $j++) { $distance += ($data[$j * $sample_number + $i] - $old_centroid[$j * $old_cluster_number + $n]) * ($data[$j * $sample_number + $i] - $old_centroid[$j * $old_cluster_number + $n]); } } $parent_variance +=$distance; } $parent_variance /=$parent_sample_number - 1; // 対数尤度を求める． $parent_likelihood =calculate_likelihood( $dimension_number, $parent_sample_number, 1, $parent_variance); // BICを求める． $parent_bic =$parent_likelihood - (($dimension_number + 1) / 2) * log($parent_sample_number); // kmeansを行う． $parent_data =array(); $counter =0; for ($i=0; $i<$sample_number; $i++) { if ($old_label[$i] == $n) { for ($j=0; $j<$dimension_number; $j++) { $parent_data[$j * $parent_sample_number + $counter] =$data[$j * $sample_number + $i]; } $counter++; } } list($child_centroid, $child_label) =kmeans_clustering( $parent_data, $dimension_number, $parent_sample_number, 2); // 新たなクラスタのデータ数を数える． $child_a_sample_number = 0; $child_b_sample_number = 0; for ($i=0; $i<$parent_sample_number; $i++) { ($child_label[$i] == 0) ? $child_a_sample_number++ : $child_b_sample_number++; } if ($child_a_sample_number <= 1 || $child_b_sample_number <= 1) { continue; } // 新たなクラスタの分散を求める． $child_a_variance =0; $child_b_variance =0; for ($i=0; $i<$parent_sample_number; $i++) { $child_a_distance = 0.0; $child_b_distance = 0.0; if ($child_label[$i] == 0) { for ($j=0; $j<$dimension_number; $j++) { $child_a_distance += ($parent_data[$j * $parent_sample_number + $i] - $child_centroid[$j * 2 + 0]) * ($parent_data[$j * $parent_sample_number + $i] - $child_centroid[$j * 2 + 0]); } } else { for ($j=0; $j<$dimension_number; $j++) { $child_b_distance += ($parent_data[$j * $parent_sample_number + $i] - $child_centroid[$j * 2 + 1]) * ($parent_data[$j * $parent_sample_number + $i] - $child_centroid[$j * 2 + 1]); } } $child_a_variance +=$child_a_distance; $child_b_variance +=$child_b_distance; } $child_a_variance /=$child_a_sample_number - 1; $child_b_variance /=$child_b_sample_number - 1; // 新たなクラスタの対数尤度を求める． $child_likelihood =calculate_likelihood( $dimension_number, $child_a_sample_number, 2, $child_a_variance ) + calculate_likelihood( $dimension_number, $child_b_sample_number, 2, $child_b_variance ); // 新たなクラスタのBICを求める． $child_bic =$child_likelihood - ((1 + 2 * $dimension_number + 2) / 2.0) * log($parent_sample_number); // 新たなクラスタのBICの方が大きければ，新たなクラスタを仲間入りさせる． if ($child_bic > $parent_bic) { // ラベルを更新する． $counter = 0; for ($i=0; $i<$sample_number; $i++) { if ($new_label[$i] == $next_label) { $new_label[$i] =$next_label + $child_label[$counter++]; } elseif ($new_label[$i] > $next_label) { $new_label[$i]++; } } // セントロイドを更新する． $temp_centroid =array_diff($new_centroid,array()); $new_centroid =array(); for ($i=0; $i<$new_cluster_number + 1; $i++) { if ($i < $next_label) { for ($j=0; $j<$dimension_number; $j++) { $new_centroid[$j * ($new_cluster_number + 1) + $i] =$temp_centroid[$j * $new_cluster_number + $i]; } } else if ($i == $next_label) { for ($j=0; $j<$dimension_number; $j++) { $new_centroid[$j * ($new_cluster_number + 1) + $i] =$child_centroid[$j * 2 + 0]; } } else if ($i == $next_label + 1) { for ($j=0; $j<$dimension_number; $j++) { $new_centroid[$j * ($new_cluster_number + 1) + $i] =$child_centroid[$j * 2 + 1]; } } else if ($i > $next_label + 1) { for ($j=0; $j<$dimension_number; $j++) { $new_centroid[$j * ($new_cluster_number + 1) + $i] =$temp_centroid[$j * $new_cluster_number + ($i - 1)]; } } } // クラスタ数を更新する． $new_cluster_number++; // 次のクラスタへ. $next_label += 2; } else { // 次のクラスタへ. $next_label++; } } $continual =$new_cluster_number != $old_cluster_number; // 結果を引き継ぐ． $old_cluster_number =$new_cluster_number; $old_centroid =array_diff($new_centroid,array()); $old_label =array_diff($new_label,array()); } while ($continual); // 最終的なクラスタを割り振る． for ($i=0; $i<$sample_number; $i++) { $label[$i] =$old_label[$i] + 1; } return array($new_centroid,$label,$new_cluster_number); } //------------------------------------- // k-means法でクラスタリング public function kmeans_clustering ( array $data, $dimension_number, $sample_number, $cluster_number) { $centroid =array(); $label =array(); $tolerance =1e-8; $old_error =PHP_INT_MAX; $new_error =PHP_INT_MAX; // ランダムにセントロイドを選択する． for ($i=0; $i<$cluster_number; $i++) { $random_id =rand(0,$sample_number-1); for ($j=0; $j<$dimension_number; $j++) { $centroid[$j * $cluster_number + $i] =$data[$j * $sample_number + $random_id]; } } // メインループ do { // 前のステップのエラー値を保持し，エラー値を初期化する． $old_error =$new_error; $new_error =0; // 前回のクラスタとセントロイドを初期化する． $cluster_size =array(); $temp_centroid =array(); for ($h=0; $h<$sample_number; $h++) { // 近いクラスタを探す． $min_distance =PHP_INT_MAX; for ($i=0; $i<$cluster_number; $i++) { $distance_from_centroid =0; for ($j=0; $j<$dimension_number; $j++) { $distance_from_centroid += ($data[$j * $sample_number + $h] - $centroid[$j * $cluster_number + $i]) * ($data[$j * $sample_number + $h] - $centroid[$j * $cluster_number + $i]); } $distance_from_centroid =sqrt($distance_from_centroid); if ($distance_from_centroid < $min_distance) { $label[$h] =$i; $min_distance =$distance_from_centroid; } } // セントロイドとクラスタ,誤差を更新する． for ($j=0; $j<$dimension_number; $j++) { $temp_centroid[$j * $cluster_number + $label[$h]] += $data[$j * $sample_number + $h]; } $cluster_size[$label[$h]]++; $new_error +=$min_distance; } // セントロイドを更新する． for ($i=0; $i<$cluster_number; $i++) { for ($j=0; $j<$dimension_number; $j++) { $centroid[$j * $cluster_number + $i] = $cluster_size[$i] ? $temp_centroid[$j * $cluster_number + $i] / $cluster_size[$i] : $temp_centroid[$j * $cluster_number + $i]; } } } while (abs($new_error - $old_error) > $tolerance); return array($centroid, $label); } //------------------------------------- // 対数尤度を求める. public function calculate_likelihood ( $dimension_number, $sample_number, $cluster_number, $variance) { $result =-(($sample_number / 2) * log(2 * pi())) -((($sample_number * $dimension_number) / 2) * log($variance)) -((($sample_number - $cluster_number) / 2)) +($sample_number * log($sample_number)); return $result; } //------------------------------------- // テスト用プロット function draw_plot ($data, $output_filename, $user_options=array()) { $options =array( "width" =>100, "height" =>100, "image_width" =>300, "image_height" =>300, "bgcolor" =>0xffffff, "plot_size" =>5 ); foreach ($user_options as $key => $value) { $options[$key] =$value; } $img =imagecreatetruecolor( $options["image_width"], $options["image_height"]); $color6hex =$options["bgcolor"]; $color =imagecolorallocate($img, ($color6hex & 0xff0000) >> 16, ($color6hex & 0x00ff00) >> 8, ($color6hex & 0x0000ff) >> 0); imagefill($img,0,0,$color); $width_factor =$options["image_width"]/$options["width"]; $height_factor =$options["image_height"]/$options["height"]; $rand_factor =rand(0,100); foreach ($data as $p => $value) { $x =$p%$options["width"]; $y =($p-$x)/$options["width"]; $color6hex =base_convert(substr(md5($value.$rand_factor),6,6),16,10); $color =imagecolorallocate($img, ($color6hex & 0xff0000) >> 16, ($color6hex & 0x00ff00) >> 8, ($color6hex & 0x0000ff) >> 0); imagefilledellipse( $img, (int)($x * $width_factor), (int)($y * $height_factor), $options["plot_size"], $options["plot_size"], $color); } ob_start(); imagepng($img); file_put_contents($output_filename,ob_get_clean()); imagedestroy($img); } /* 1辺100の正方形(2次元)内で12点サンプル。 初期クラスタ数5。 */ $size =100; $cluster_number =5; $sample_number =12; $dimension =2; $data =array(); $data[0*$dimension+0] =9; $data[1*$dimension+0] =9; $data[0*$dimension+1] =15; $data[1*$dimension+1] =11; $data[0*$dimension+2] =11; $data[1*$dimension+2] =14; $data[0*$dimension+3] =12; $data[1*$dimension+3] =13; $data[0*$dimension+4] =43; $data[1*$dimension+4] =43; $data[0*$dimension+5] =48; $data[1*$dimension+5] =44; $data[0*$dimension+6] =86; $data[1*$dimension+6] =57; $data[0*$dimension+7] =87; $data[1*$dimension+7] =48; $data[0*$dimension+8] =36; $data[1*$dimension+8] =37; $data[0*$dimension+9] =37; $data[1*$dimension+9] =30; $data[0*$dimension+10] =55; $data[1*$dimension+10] =48; $data[0*$dimension+11] =56; $data[1*$dimension+11] =40; list($centroid,$label,$new_cluster_number) =xmeans_clustering($data,$dimension,$sample_number,$cluster_number); $plot =array(); for ($i=0; $i<$sample_number; $i++) { $x =$data[0*$dimension+$i]; $y =$data[1*$dimension+$i]; $plot[$y*$size+$x] =$label[$i]; } draw_plot($plot,"./dest.png",array( "width" =>$size, "height" =>$size, "bgcolor" =>0xeeeeee, "plot_size" =>10, )); print '<html><body><img src="./dest.png"></body></html>'; ?>
コードばっかりですが、ご参考まで。
thx
