クラスター分析

/****************************/
/* クラスター分析           */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>

void cluster(int, int, int, int, double **, int *);
double range(int, int, int, int, double **);
double range_c(int, double, double, double, int, int, int);

int main()
{
	double **x;
	int i1, i2, k, k1, L, n, N, *g;
	int method;
					// 方法,グループ数,変数の数,及び,データの数
	scanf("%d %d %d %d", &method, &L, &n, &N);

	g = new int [N];
	x = new double * [N];
					// データ
	for (i1 = 0; i1 < N; i1++) {
		x[i1] = new double [n];
		for (i2 = 0; i2 < n; i2++)
			scanf("%lf", &x[i1][i2]);
	}
					// クラスター分析
	cluster(method, L, N, n, x, g);
					// 出力
	k = 1;
	while (k <= L) {
		k1 = -1;
		printf("グループ %d\n", k);
		for (i1 = 0; i1 < N && k1 < 0; i1++) {
			if (g[i1] >= 0)
				k1 = g[i1];
		}
		for (i1 = 0; i1 < N; i1++) {
			if (g[i1] == k1) {
				printf("   %d", i1);
				for (i2 = 0; i2 < n; i2++)
					printf(" %f", x[i1][i2]);
				printf("\n");
				g[i1] = -1;
			}
		}
		k++;
	}

	for (i1 = 0; i1 < N; i1++)
		delete [] x[i1];
	delete [] x;
	delete [] g;

	return 0;
}

/*********************************/
/* クラスター分析                */
/*      method : =1 : 最短距離法 */
/*               =2 : 最長距離法 */
/*               =3 : メジアン法 */
/*               =4 : 重心法     */
/*               =5 : 群平均法   */
/*               =6 : ウォード法 */
/*      L : グループの数         */
/*      N : データの数           */
/*      n : 変量の数             */
/*      x : データ               */
/*      g : 所属するグループ番号 */
/*********************************/
#include <math.h>

void cluster(int method, int L, int N, int n, double **x, int *g)
{
	double **r, min;
	int ci = 0, cj = 0, i1, i2, k, M = N, *n_g;
					// 初期設定
	k   = (method < 4) ? 0 : 1;
	n_g = new int [N];
	r   = new double * [N];
	for (i1 = 0; i1 < N; i1++) {
		g[i1]   = i1;
		n_g[i1] = 1;
		r[i1]   = new double [N];
		for (i2 = i1+1; i2 < N; i2++)
			r[i1][i2] = range(k, i1, i2, n, x);
	}
	for (i1 = 0; i1 < N; i1++) {
		for (i2 = i1+1; i2 < N; i2++)
			r[i2][i1] = r[i1][i2];
	}
					// 実行
	while (M > L) {
						// 最小距離のクラスターを探す
		min = -1.0;
		for (i1 = 0; i1 < N; i1++) {
			if (g[i1] == i1) {
				for (i2 = i1+1; i2 < N; i2++) {
					if (g[i2] == i2) {
						if (min < 0.0 || r[i1][i2] < min) {
							min = r[i1][i2];
							ci  = i1;
							cj  = i2;
						}
					}
				}
			}
		}
						// クラスターを融合し,他のクラスターとの距離を計算
		for (i1 = 0; i1 < N; i1++) {
			if (g[i1] == i1) {
				for (i2 = i1+1; i2 < N; i2++) {
					if (g[i2] == i2) {
						if (i1 != cj && i2 != cj) {
							if (i1 == ci) {
								r[i1][i2] = range_c(method, r[ci][i2], r[cj][i2], r[ci][cj],
                                                    n_g[ci], n_g[cj], n_g[i2]);
								r[i2][i1] = r[i1][i2];
							}
							else if (i2 == ci) {
								r[i1][i2] = range_c(method, r[i1][ci], r[i1][cj], r[ci][cj],
                                                    n_g[ci], n_g[cj], n_g[i2]);
								r[i2][i1] = r[i1][i2];
							}
						}
					}
				}
			}
		}

		n_g[ci] += n_g[cj];
		for (i1 = 0; i1 < N; i1++) {
			if (g[i1] == cj)
				g[i1] = ci;
		}

		M--;
	}
					// 領域の開放
	for (i1 = 0; i1 < N; i1++)
		delete [] r[i1];
	delete [] r;
	delete [] n_g;
}

/*********************************************/
/* 2つのデータ間の距離                      */
/*      method : =0 : ユークリッド距離       */
/*               =1 : ユークリッド距離の2乗 */
/*      i,j : データ番号                     */
/*      n : 変量の数                         */
/*      x : データ                           */
/*      return : 距離                        */
/*********************************************/
double range(int method, int i, int j, int n, double **x)
{
	double x1, r = 0.0;
	int i2;

	for (i2 = 0; i2 < n; i2++) {
		x1  = x[i][i2] - x[j][i2];
		r  += x1 * x1;
	}
	if (method == 0)
		r = sqrt(r);

	return r;
}

/***********************************************/
/* 2つのクラスター間の距離                    */
/*      method : =1 : 最短距離法               */
/*               =2 : 最長距離法               */
/*               =3 : メジアン法               */
/*               =4 : 重心法                   */
/*               =5 : 群平均法                 */
/*               =6 : ウォード法               */
/*      Dio : クラスターiとクラスターoとの距離 */
/*      Djo : クラスターjとクラスターoとの距離 */
/*      Dij : クラスターiとクラスターjとの距離 */
/*      ni : クラスターiに含まれるデータ数     */
/*      nj : クラスターjに含まれるデータ数     */
/*      no : クラスターoに含まれるデータ数     */
/*      x,y : データ                           */
/*      return : 距離                          */
/***********************************************/
double range_c(int method, double Dio, double Djo, double Dij, int ni, int nj, int no)
{
	double r = 0.0;
	int nk;

	switch (method) {
					// 最短距離法
		case 1:
			r = (Dio <= Djo) ? Dio : Djo;
			break;
					// 最長距離法
		case 2:
			r = (Dio >= Djo) ? Dio : Djo;
			break;
					// メジアン法
		case 3:
			r = 0.5 * Dio + 0.5 * Djo - 0.25 * Dij;
			break;
					// 重心法
		case 4:
			nk = ni + nj;
			r  = ni * Dio / nk + nj * Djo / nk - (double)ni * nj * Dij / ((double)nk * nk);
			break;
					// 群平均法
		case 5:
			nk = ni + nj;
			r  = ni * Dio / nk + nj * Djo / nk;
			break;
					// ウォード法
		case 6:
			nk = ni + nj;
			r  = ((ni + no) * Dio + (nj + no) * Djo - no * Dij) / ((double)nk + no);
			break;
	}

	return r;
}

---------データ例(コメント部分を除いて下さい)---------
1 3 2 100   // クラスター間距離を計算する方法(method),クラスタの数(L,クラスタの数がこの値になったら終了),変数の数(n),及び,データの数(N)
66.813742 25.509139   // x, y
50.813965 25.537399
74.048072 39.627578
14.792737 67.086620
42.838304 11.391394
32.301288 68.001189
24.436877 42.170217
60.451432 36.760112
24.314984 30.052427
55.573736 63.561140
72.162565 19.845151
43.217518 90.320780
51.888313 53.074269
28.150842 40.674284
31.619668 57.984410
64.360526 49.975656
53.702823 79.538952
45.052136 35.297650
4.702717 49.420025
92.649989 5.431771
65.403044 56.028507
48.797150 30.888271
40.806244 14.484027
84.355436 53.896156
50.133508 76.624964
49.405414 78.403179
76.824779 41.008448
63.246400 40.075843
78.538747 62.910356
64.396460 56.299403
77.152853 44.325319
80.592411 36.177862
48.763290 43.966714
58.789819 90.533976
13.684999 67.266671
75.163001 81.909246
46.451671 56.376068
65.105023 34.563928
45.815747 20.164535
33.227702 44.477025
41.061314 -54.334662
40.021202 -65.490282
22.829721 -56.861792
43.917033 -43.558827
40.723550 -43.188580
53.309317 -69.996428
31.463381 -31.506555
24.313168 -52.662580
17.433859 -57.468181
32.457761 -93.127719
9.183688 -31.686863
10.186332 -44.173216
27.925638 -62.236940
9.356182 -42.291381
23.036540 -27.153692
-16.261328 -35.185462
29.557288 -31.933335
36.519228 -59.269642
25.692837 -57.033021
22.797337 -52.169044
15.312205 -70.370841
22.086548 -25.219727
50.448475 -34.825824
30.814475 -40.968689
28.798789 -65.737936
59.822791 -25.673863
16.224354 -48.408079
6.867234 -61.804123
44.250439 -50.598095
54.117605 -20.539278
-43.244836 -3.374350
-19.122981 -9.026765
-4.632598 27.480272
26.647517 24.535738
19.527239 1.573533
-29.982384 26.884977
-11.856633 -24.568687
-9.244616 -3.343794
-15.656771 32.298241
-22.152364 24.203428
-34.526986 -27.862305
-60.500060 -18.400603
-6.902336 21.631805
-46.863390 -27.364383
-36.999128 4.181449
-24.599285 -17.177658
-24.814923 9.906257
-7.845930 9.388576
-36.128527 5.644606
-31.522759 11.282073
-12.500001 24.953627
-21.330455 -0.830663
-21.270652 -1.890803
-23.750951 5.435622
-25.563039 22.105947
-16.047432 -5.815069
-34.280972 30.657254
-36.264330 -8.759981
-12.109507 -15.804752
-33.227330 -5.857611