静岡理工科大学 菅沼ホーム 目次 索引

主成分分析

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

  プログラムは,主成分分析を行うためのものです.データの入力方法に関しては,プログラムの最後に提示してある入力例に対する説明を参考にして下さい.

  1. C++

    /****************************/
    /* 主成分分析               */
    /*      coded by Y.Suganuma */
    /****************************/
    #include <stdio.h>
    #include <math.h>
    
    int principal(int, int, double **, double *, double **, double, int);
    int Jacobi(int, int, double, double **, double **, double **, double **, double **);
    
    int main()
    {
    	double **x, *r, **a;
    	int i1, i2, n, p, sw;
    
    	scanf("%d %d", &p, &n);   // 変数の数とデータの数
    
    	r = new double [p];
    	x = new double * [p];
    	a = new double * [p];
    	for (i1 = 0; i1 < p; i1++) {
    		x[i1]  = new double [n];
    		a[i1]  = new double [p];
    	}
    
    	for (i1 = 0; i1 < n; i1++) {   // データ
    		for (i2 = 0; i2 < p; i2++)
    			scanf("%lf", &x[i2][i1]);
    	}
    
    	sw = principal(p, n, x, r, a, 1.0e-10, 200);
    
    	if (sw == 0) {
    		for (i1 = 0; i1 < p; i1++) {
    			printf("主成分 %f", r[i1]);
    			printf(" 係数");
    			for (i2 = 0; i2 < p; i2++)
    				printf(" %f", a[i1][i2]);
    			printf("\n");
    		}
    	}
    	else
    		printf("***error  解を求めることができませんでした\n");
    
    	for (i1 = 0; i1 < p; i1++) {
    		delete [] x[i1];
    		delete [] a[i1];
    	}
    	delete [] x;
    	delete [] a;
    	delete [] r;
    
    	return 0;
    }
    
    /***********************************/
    /* 主成分分析                      */
    /*      p : 変数の数               */
    /*      n : データの数             */
    /*      x : データ                 */
    /*      r : 分散(主成分)         */
    /*      a : 係数                   */
    /*      eps : 収束性を判定する規準 */
    /*      ct : 最大繰り返し回数      */
    /*      return : =0 : 正常         */
    /*               =1 : エラー       */
    /***********************************/
    int principal(int p, int n, double **x, double *r, double **a, double eps, int ct)
    {
    	double **A1, **A2, **C, mean, **X1, **X2, s2;
    	int i1, i2, i3, sw = 0;
    					// 領域の確保
    	C  = new double * [p];
    	A1 = new double * [p];
    	A2 = new double * [p];
    	X1 = new double * [p];
    	X2 = new double * [p];
    	for (i1 = 0; i1 < p; i1++) {
    		C[i1]  = new double [p];
    		A1[i1] = new double [p];
    		A2[i1] = new double [p];
    		X1[i1] = new double [p];
    		X2[i1] = new double [p];
    	}
    					// データの基準化
    	for (i1 = 0; i1 < p; i1++) {
    		mean = 0.0;
    		s2   = 0.0;
    		for (i2 = 0; i2 < n; i2++) {
    			mean += x[i1][i2];
    			s2   += x[i1][i2] * x[i1][i2];
    		}
    		mean /= n;
    		s2   /= n;
    		s2    = n * (s2 - mean * mean) / (n - 1);
    		s2    = sqrt(s2);
    		for (i2 = 0; i2 < n; i2++)
    			x[i1][i2] = (x[i1][i2] - mean) / s2;
    	}
    					// 分散強分散行列の計算
    	for (i1 = 0; i1 < p; i1++) {
    		for (i2 = i1; i2 < p; i2++) {
    			s2 = 0.0;
    			for (i3 = 0; i3 < n; i3++)
    				s2 += x[i1][i3] * x[i2][i3];
    			s2        /= (n - 1);
    			C[i1][i2]  = s2;
    			if (i1 != i2)
    				C[i2][i1] = s2;
    		}
    	}
    					// 固有値と固有ベクトルの計算(ヤコビ法)
    	sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
    
    	if (sw == 0) {
    		for (i1 = 0; i1 < p; i1++) {
    			r[i1] = A1[i1][i1];
    			for (i2 = 0; i2 < p; i2++)
    				a[i1][i2] = X1[i2][i1];
    		}
    	}
    					// 領域の解放
    	for (i1 = 0; i1 < p; i1++) {
    		delete [] C[i1];
    		delete [] A1[i1];
    		delete [] A2[i1];
    		delete [] X1[i1];
    		delete [] X2[i1];
    	}
    	delete [] C;
    	delete [] A1;
    	delete [] A2;
    	delete [] X1;
    	delete [] X2;
    
    	return sw;
    }
    
    /*************************************************************/
    /* 実対称行列の固有値・固有ベクトル(ヤコビ法)              */
    /*      n : 次数                                             */
    /*      ct : 最大繰り返し回数                                */
    /*      eps : 収束判定条件                                   */
    /*      A : 対象とする行列                                   */
    /*      A1, A2 : 作業域(nxnの行列),A1の対角要素が固有値   */
    /*      X1, X2 : 作業域(nxnの行列),X1の各列が固有ベクトル */
    /*      return : =0 : 正常                                   */
    /*               =1 : 収束せず                               */
    /*      coded by Y.Suganuma                                  */
    /*************************************************************/
    #include <math.h>
    
    int Jacobi(int n, int ct, double eps, double **A, double **A1, double **A2,
               double **X1, double **X2)
    {
    	double max, s, t, v, sn, cs;
    	int i1, i2, k = 0, ind = 1, p = 0, q = 0;
    					// 初期設定
    	for (i1 = 0; i1 < n; i1++) {
    		for (i2 = 0; i2 < n; i2++) {
    			A1[i1][i2] = A[i1][i2];
    			X1[i1][i2] = 0.0;
    		}
    		X1[i1][i1] = 1.0;
    	}
    					// 計算
    	while (ind > 0 && k < ct) {
    						// 最大要素の探索
    		max = 0.0;
    		for (i1 = 0; i1 < n; i1++) {
    			for (i2 = 0; i2 < n; i2++) {
    				if (i2 != i1) {
    					if (fabs(A1[i1][i2]) > max) {
    						max = fabs(A1[i1][i2]);
    						p   = i1;
    						q   = i2;
    					}
    				}
    			}
    		}
    						// 収束判定
    							// 収束した
    		if (max < eps)
    			ind = 0;
    							// 収束しない
    		else {
    								// 準備
    			s  = -A1[p][q];
    			t  = 0.5 * (A1[p][p] - A1[q][q]);
    			v  = fabs(t) / sqrt(s * s + t * t);
    			sn = sqrt(0.5 * (1.0 - v));
    			if (s*t < 0.0)
    				sn = -sn;
    			cs = sqrt(1.0 - sn * sn);
    								// Akの計算
    			for (i1 = 0; i1 < n; i1++) {
    				if (i1 == p) {
    					for (i2 = 0; i2 < n; i2++) {
    						if (i2 == p)
    							A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                           2.0 * A1[p][q] * sn * cs;
    						else if (i2 == q)
    							A2[p][q] = 0.0;
    						else
    							A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
    					}
    				}
    				else if (i1 == q) {
    					for (i2 = 0; i2 < n; i2++) {
    						if (i2 == q)
    							A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                           2.0 * A1[p][q] * sn * cs;
    						else if (i2 == p)
    							A2[q][p] = 0.0;
    						else
    							A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
    					}
    				}
    				else {
    					for (i2 = 0; i2 < n; i2++) {
    						if (i2 == p)
    							A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
    						else if (i2 == q)
    							A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
    						else
    							A2[i1][i2] = A1[i1][i2];
    					}
    				}
    			}
    								// Xkの計算
    			for (i1 = 0; i1 < n; i1++) {
    				for (i2 = 0; i2 < n; i2++) {
    					if (i2 == p)
    						X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
    					else if (i2 == q)
    						X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
    					else
    						X2[i1][i2] = X1[i1][i2];
    				}
    			}
    								// 次のステップへ
    			k++;
    			for (i1 = 0; i1 < n; i1++) {
    				for (i2 = 0; i2 < n; i2++) {
    					A1[i1][i2] = A2[i1][i2];
    					X1[i1][i2] = X2[i1][i2];
    				}
    			}
    		}
    	}
    
    	return ind;
    }
    
    /*
    ---------データ例(コメント部分を除いて下さい)---------
    4 100   // 変数の数(p)とデータの数(n)
    66 22 44 31   // x1, x2, x3, x4
    25 74 17 81
    50 23 53 71
    25 57 19 81
    74 47 64 47
    39 33 48 46
    14 22 9 69
    67 60 49 26
    42 40 77 65
    11 80 0 86
    32 0 43 74
    68 69 44 68
    24 49 9 71
    42 74 28 46
    60 58 73 28
    36 37 33 68
    24 44 19 83
    30 40 31 50
    55 40 60 49
    63 47 94 41
    72 30 100 45
    19 22 13 75
    43 39 43 34
    90 83 92 31
    51 77 52 82
    53 70 34 31
    28 51 53 44
    40 62 42 79
    31 48 22 68
    57 29 51 30
    64 89 57 42
    49 82 72 29
    53 31 55 43
    79 52 70 10
    45 19 43 57
    35 34 34 89
    4 69 0 100
    49 49 66 66
    92 82 97 6
    5 89 0 100
    65 26 83 28
    56 36 64 38
    48 50 25 22
    30 30 15 55
    40 65 38 42
    14 67 9 67
    84 96 90 8
    53 64 51 54
    50 89 60 52
    76 41 68 9
    49 40 53 49
    78 66 66 17
    76 58 90 29
    41 15 40 49
    63 60 55 33
    40 36 49 67
    78 54 71 18
    62 72 69 12
    64 47 42 53
    56 64 9 15
    77 35 56 25
    44 12 46 87
    80 9 56 19
    36 21 52 78
    48 63 64 48
    43 61 50 47
    58 23 28 50
    90 12 100 0
    13 33 11 77
    67 44 48 28
    75 45 68 17
    81 22 89 9
    46 45 59 55
    56 49 64 55
    65 62 72 27
    34 49 29 77
    45 33 60 63
    20 45 14 99
    33 38 26 87
    44 51 69 52
    64 57 64 48
    44 64 51 28
    63 48 56 11
    29 39 33 84
    40 48 51 54
    40 38 26 62
    68 46 61 26
    58 45 68 48
    64 44 77 63
    59 62 44 66
    81 53 93 19
    23 34 12 68
    51 35 55 46
    74 70 84 17
    42 33 56 44
    46 31 46 53
    33 57 38 63
    40 24 20 42
    53 36 60 31
    0 34 0 100
    */
    			

  2. Java

      アプレット版では,任意のデータに対して画面上で実行することができます.
    /****************************/
    /* 主成分分析               */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.io.*;
    import java.util.StringTokenizer;
    
    public class Test {
    	public static void main(String args[]) throws IOException
    	{
    		double x[][], r[], a[][];
    		int i1, i2, n, p, sw;
    		StringTokenizer str;
    		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
    					// 変数の数とデータの数
    		str = new StringTokenizer(in.readLine(), " ");
    		p   = Integer.parseInt(str.nextToken());
    		n   = Integer.parseInt(str.nextToken());
    
    		r = new double [p];
    		x = new double [p][n];
    		a = new double [p][p];
    					// データ
    		for (i1 = 0; i1 < n; i1++) {
    			str = new StringTokenizer(in.readLine(), " ");
    			for (i2 = 0; i2 < p; i2++)
    				x[i2][i1] = Double.parseDouble(str.nextToken());
    		}
    
    		sw = principal(p, n, x, r, a, 1.0e-10, 200);
    
    		if (sw == 0) {
    			for (i1 = 0; i1 < p; i1++) {
    				System.out.print("主成分 " + r[i1]);
    				System.out.print(" 係数");
    				for (i2 = 0; i2 < p; i2++)
    					System.out.print(" " + a[i1][i2]);
    				System.out.println();
    			}
    		}
    		else
    			System.out.println("***error  解を求めることができませんでした");
    	}
    
    	/***********************************/
    	/* 主成分分析                      */
    	/*      p : 変数の数               */
    	/*      n : データの数             */
    	/*      x : データ                 */
    	/*      r : 分散(主成分)         */
    	/*      a : 係数                   */
    	/*      eps : 収束性を判定する規準 */
    	/*      ct : 最大繰り返し回数      */
    	/*      return : =0 : 正常         */
    	/*               =1 : エラー       */
    	/***********************************/
    	static int principal(int p, int n, double x[][], double r[], double a[][], double eps, int ct)
    	{
    		double A1[][], A2[][], C[][], mean, X1[][], X2[][], s2;
    		int i1, i2, i3, sw = 0;
    					// 領域の確保
    		C  = new double [p][p];
    		A1 = new double [p][p];
    		A2 = new double [p][p];
    		X1 = new double [p][p];
    		X2 = new double [p][p];
    					// データの基準化
    		for (i1 = 0; i1 < p; i1++) {
    			mean = 0.0;
    			s2   = 0.0;
    			for (i2 = 0; i2 < n; i2++) {
    				mean += x[i1][i2];
    				s2   += x[i1][i2] * x[i1][i2];
    			}
    			mean /= n;
    			s2   /= n;
    			s2    = n * (s2 - mean * mean) / (n - 1);
    			s2    = Math.sqrt(s2);
    			for (i2 = 0; i2 < n; i2++)
    				x[i1][i2] = (x[i1][i2] - mean) / s2;
    		}
    					// 分散強分散行列の計算
    		for (i1 = 0; i1 < p; i1++) {
    			for (i2 = i1; i2 < p; i2++) {
    				s2 = 0.0;
    				for (i3 = 0; i3 < n; i3++)
    					s2 += x[i1][i3] * x[i2][i3];
    				s2        /= (n - 1);
    				C[i1][i2]  = s2;
    				if (i1 != i2)
    					C[i2][i1] = s2;
    			}
    		}
    					// 固有値と固有ベクトルの計算(ヤコビ法)
    		sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
    
    		if (sw == 0) {
    			for (i1 = 0; i1 < p; i1++) {
    				r[i1] = A1[i1][i1];
    				for (i2 = 0; i2 < p; i2++)
    					a[i1][i2] = X1[i2][i1];
    			}
    		}
    
    		return sw;
    	}
    
    	/*************************************************************/
    	/* 実対称行列の固有値・固有ベクトル(ヤコビ法)              */
    	/*      n : 次数                                             */
    	/*      ct : 最大繰り返し回数                                */
    	/*      eps : 収束判定条件                                   */
    	/*      A : 対象とする行列                                   */
    	/*      A1, A2 : 作業域(nxnの行列),A1の対角要素が固有値   */
    	/*      X1, X2 : 作業域(nxnの行列),X1の各列が固有ベクトル */
    	/*      return : =0 : 正常                                   */
    	/*               =1 : 収束せず                               */
    	/*      coded by Y.Suganuma                                  */
    	/*************************************************************/
    	static int Jacobi(int n, int ct, double eps, double A[][], double A1[][], double A2[][],
    	                  double X1[][], double X2[][])
    	{
    		double max, s, t, v, sn, cs;
    		int i1, i2, k = 0, ind = 1, p = 0, q = 0;
    					// 初期設定
    		for (i1 = 0; i1 < n; i1++) {
    			for (i2 = 0; i2 < n; i2++) {
    				A1[i1][i2] = A[i1][i2];
    				X1[i1][i2] = 0.0;
    			}
    			X1[i1][i1] = 1.0;
    		}
    					// 計算
    		while (ind > 0 && k < ct) {
    						// 最大要素の探索
    			max = 0.0;
    			for (i1 = 0; i1 < n; i1++) {
    				for (i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						if (Math.abs(A1[i1][i2]) > max) {
    							max = Math.abs(A1[i1][i2]);
    							p   = i1;
    							q   = i2;
    						}
    					}
    				}
    			}
    						// 収束判定
    							// 収束した
    			if (max < eps)
    				ind = 0;
    							// 収束しない
    			else {
    								// 準備
    				s  = -A1[p][q];
    				t  = 0.5 * (A1[p][p] - A1[q][q]);
    				v  = Math.abs(t) / Math.sqrt(s * s + t * t);
    				sn = Math.sqrt(0.5 * (1.0 - v));
    				if (s*t < 0.0)
    					sn = -sn;
    				cs = Math.sqrt(1.0 - sn * sn);
    								// Akの計算
    				for (i1 = 0; i1 < n; i1++) {
    					if (i1 == p) {
    						for (i2 = 0; i2 < n; i2++) {
    							if (i2 == p)
    								A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                               2.0 * A1[p][q] * sn * cs;
    							else if (i2 == q)
    								A2[p][q] = 0.0;
    							else
    								A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
    						}
    					}
    					else if (i1 == q) {
    						for (i2 = 0; i2 < n; i2++) {
    							if (i2 == q)
    								A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                               2.0 * A1[p][q] * sn * cs;
    							else if (i2 == p)
    								A2[q][p] = 0.0;
    							else
    								A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
    						}
    					}
    					else {
    						for (i2 = 0; i2 < n; i2++) {
    							if (i2 == p)
    								A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
    							else if (i2 == q)
    								A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
    							else
    								A2[i1][i2] = A1[i1][i2];
    						}
    					}
    				}
    								// Xkの計算
    				for (i1 = 0; i1 < n; i1++) {
    					for (i2 = 0; i2 < n; i2++) {
    						if (i2 == p)
    							X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
    						else if (i2 == q)
    							X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
    						else
    							X2[i1][i2] = X1[i1][i2];
    					}
    				}
    								// 次のステップへ
    				k++;
    				for (i1 = 0; i1 < n; i1++) {
    					for (i2 = 0; i2 < n; i2++) {
    						A1[i1][i2] = A2[i1][i2];
    						X1[i1][i2] = X2[i1][i2];
    					}
    				}
    			}
    		}
    
    		return ind;
    	}
    }
    
    /*
    ---------データ例(コメント部分を除いて下さい)---------
    4 100   // 変数の数(p)とデータの数(n)
    66 22 44 31   // x1, x2, x3, x4
    25 74 17 81
    50 23 53 71
    25 57 19 81
    74 47 64 47
    39 33 48 46
    14 22 9 69
    67 60 49 26
    42 40 77 65
    11 80 0 86
    32 0 43 74
    68 69 44 68
    24 49 9 71
    42 74 28 46
    60 58 73 28
    36 37 33 68
    24 44 19 83
    30 40 31 50
    55 40 60 49
    63 47 94 41
    72 30 100 45
    19 22 13 75
    43 39 43 34
    90 83 92 31
    51 77 52 82
    53 70 34 31
    28 51 53 44
    40 62 42 79
    31 48 22 68
    57 29 51 30
    64 89 57 42
    49 82 72 29
    53 31 55 43
    79 52 70 10
    45 19 43 57
    35 34 34 89
    4 69 0 100
    49 49 66 66
    92 82 97 6
    5 89 0 100
    65 26 83 28
    56 36 64 38
    48 50 25 22
    30 30 15 55
    40 65 38 42
    14 67 9 67
    84 96 90 8
    53 64 51 54
    50 89 60 52
    76 41 68 9
    49 40 53 49
    78 66 66 17
    76 58 90 29
    41 15 40 49
    63 60 55 33
    40 36 49 67
    78 54 71 18
    62 72 69 12
    64 47 42 53
    56 64 9 15
    77 35 56 25
    44 12 46 87
    80 9 56 19
    36 21 52 78
    48 63 64 48
    43 61 50 47
    58 23 28 50
    90 12 100 0
    13 33 11 77
    67 44 48 28
    75 45 68 17
    81 22 89 9
    46 45 59 55
    56 49 64 55
    65 62 72 27
    34 49 29 77
    45 33 60 63
    20 45 14 99
    33 38 26 87
    44 51 69 52
    64 57 64 48
    44 64 51 28
    63 48 56 11
    29 39 33 84
    40 48 51 54
    40 38 26 62
    68 46 61 26
    58 45 68 48
    64 44 77 63
    59 62 44 66
    81 53 93 19
    23 34 12 68
    51 35 55 46
    74 70 84 17
    42 33 56 44
    46 31 46 53
    33 57 38 63
    40 24 20 42
    53 36 60 31
    0 34 0 100
    */
    			
    アプレット版
    /****************************/
    /* 主成分分析               */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.awt.*;
    import java.awt.event.*;
    import java.applet.*;
    import java.util.StringTokenizer;
    
    public class Prin extends Applet implements ActionListener {
    	TextField tx1, tx2;
    	TextArea ta1, ta2;
    	Button bt;
    	Font f = new Font("TimesRoman", Font.BOLD, 20);
    
    	/************/
    	/* 初期設定 */
    	/************/
    	public void init()
    	{
    					// レイアウト,背景色,フォント
    		setLayout(new BorderLayout(5, 5));
    		setBackground(new Color(225, 255, 225));
    		setFont(f);
    					// 上のパネル
    		Panel pn1 = new Panel();
    		add(pn1, BorderLayout.NORTH);
    
    		pn1.add(new Label("変数の数:"));
    		tx1 = new TextField("4", 2);
    		pn1.add(tx1);
    
    		pn1.add(new Label("データの数:"));
    		tx2 = new TextField("100", 3);
    		pn1.add(tx2);
    					// 中央のパネル(データ)
    		Panel pn2 = new Panel();
    		add(pn2, BorderLayout.CENTER);
    
    		pn2.add(new Label("データ:"));
    		ta1 = new TextArea("66 22 44 31\n25 74 17 81\n50 23 53 71\n25 57 19 81\n74 47 64 47\n39 33 48 46\n14 22 9 69\n67 60 49 26\n42 40 77 65\n11 80 0 86\n32 0 43 74\n68 69 44 68\n24 49 9 71\n42 74 28 46\n60 58 73 28\n36 37 33 68\n24 44 19 83\n30 40 31 50\n55 40 60 49\n63 47 94 41\n72 30 100 45\n19 22 13 75\n43 39 43 34\n90 83 92 31\n51 77 52 82\n53 70 34 31\n28 51 53 44\n40 62 42 79\n31 48 22 68\n57 29 51 30\n64 89 57 42\n49 82 72 29\n53 31 55 43\n79 52 70 10\n45 19 43 57\n35 34 34 89\n4 69 0 100\n49 49 66 66\n92 82 97 6\n5 89 0 100\n65 26 83 28\n56 36 64 38\n48 50 25 22\n30 30 15 55\n40 65 38 42\n14 67 9 67\n84 96 90 8\n53 64 51 54\n50 89 60 52\n76 41 68 9\n49 40 53 49\n78 66 66 17\n76 58 90 29\n41 15 40 49\n63 60 55 33\n40 36 49 67\n78 54 71 18\n62 72 69 12\n64 47 42 53\n56 64 9 15\n77 35 56 25\n44 12 46 87\n80 9 56 19\n36 21 52 78\n48 63 64 48\n43 61 50 47\n58 23 28 50\n90 12 100 0\n13 33 11 77\n67 44 48 28\n75 45 68 17\n81 22 89 9\n46 45 59 55\n56 49 64 55\n65 62 72 27\n34 49 29 77\n45 33 60 63\n20 45 14 99\n33 38 26 87\n44 51 69 52\n64 57 64 48\n44 64 51 28\n63 48 56 11\n29 39 33 84\n40 48 51 54\n40 38 26 62\n68 46 61 26\n58 45 68 48\n64 44 77 63\n59 62 44 66\n81 53 93 19\n23 34 12 68\n51 35 55 46\n74 70 84 17\n42 33 56 44\n46 31 46 53\n33 57 38 63\n40 24 20 42\n53 36 60 31\n0 34 0 100\n", 15, 15);
    		pn2.add(ta1);
    
    		pn2.add(new Label(" "));
    		bt = new Button("OK");
    		bt.setBackground(Color.pink);
    		bt.addActionListener(this);
    		pn2.add(bt);
    					// 下のパネル(結果)
    		Panel pn3 = new Panel();
    		add(pn3, BorderLayout.SOUTH);
    
    		pn3.add(new Label("結果:"));
    		ta2 = new TextArea(7, 70);
    		pn3.add(ta2);
    	}
    
    	/******************************/
    	/* ボタンが押されたときの処理 */
    	/******************************/
    	public void actionPerformed(ActionEvent e)
    	{
    		if (e.getSource() == bt) {
    					// データ
    			int p = Integer.parseInt(tx1.getText());
    			int n = Integer.parseInt(tx2.getText());
    
    			double r[] = new double [p];
    			double x[][] = new double [p][n];
    			double a[][] = new double [p][p];
    
    			String ss = ta1.getText();
    			StringTokenizer str = new StringTokenizer(ss, " \n");
    			for (int i1 = 0; str.hasMoreTokens() && i1 < n; i1++) {
    				for (int i2 = 0; i2 < p; i2++)
    					x[i2][i1] = Double.parseDouble(str.nextToken());
    			}
    					// 計算
    			int sw = principal(p, n, x, r, a, 1.0e-10, 200);
    
    			if (sw > 0)
    				ta2.setText("収束しません");
    			else {
    				ta2.setText("");
    				for (int i1 = 0; i1 < p; i1++) {
    					ta2.append("主成分 " + String.format("%.5f",r[i1]));
    					ta2.append(" 係数");
    					for (int i2 = 0; i2 < p; i2++)
    						ta2.append(" " + String.format("%.5f",a[i1][i2]));
    					ta2.append("\n");
    				}
    			}
    		}
    	}
    
    	/***********************************/
    	/* 主成分分析                      */
    	/*      p : 変数の数               */
    	/*      n : データの数             */
    	/*      x : データ                 */
    	/*      r : 分散(主成分)         */
    	/*      a : 係数                   */
    	/*      eps : 収束性を判定する規準 */
    	/*      ct : 最大繰り返し回数      */
    	/*      return : =0 : 正常         */
    	/*               =1 : エラー       */
    	/***********************************/
    	int principal(int p, int n, double x[][], double r[], double a[][], double eps, int ct)
    	{
    		double A1[][], A2[][], C[][], mean, X1[][], X2[][], s2;
    		int i1, i2, i3, sw = 0;
    					// 領域の確保
    		C  = new double [p][p];
    		A1 = new double [p][p];
    		A2 = new double [p][p];
    		X1 = new double [p][p];
    		X2 = new double [p][p];
    					// データの基準化
    		for (i1 = 0; i1 < p; i1++) {
    			mean = 0.0;
    			s2   = 0.0;
    			for (i2 = 0; i2 < n; i2++) {
    				mean += x[i1][i2];
    				s2   += x[i1][i2] * x[i1][i2];
    			}
    			mean /= n;
    			s2   /= n;
    			s2    = n * (s2 - mean * mean) / (n - 1);
    			s2    = Math.sqrt(s2);
    			for (i2 = 0; i2 < n; i2++)
    				x[i1][i2] = (x[i1][i2] - mean) / s2;
    		}
    					// 分散強分散行列の計算
    		for (i1 = 0; i1 < p; i1++) {
    			for (i2 = i1; i2 < p; i2++) {
    				s2 = 0.0;
    				for (i3 = 0; i3 < n; i3++)
    					s2 += x[i1][i3] * x[i2][i3];
    				s2        /= (n - 1);
    				C[i1][i2]  = s2;
    				if (i1 != i2)
    					C[i2][i1] = s2;
    			}
    		}
    					// 固有値と固有ベクトルの計算(ヤコビ法)
    		sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
    
    		if (sw == 0) {
    			for (i1 = 0; i1 < p; i1++) {
    				r[i1] = A1[i1][i1];
    				for (i2 = 0; i2 < p; i2++)
    					a[i1][i2] = X1[i2][i1];
    			}
    		}
    
    		return sw;
    	}
    
    	/*************************************************************/
    	/* 実対称行列の固有値・固有ベクトル(ヤコビ法)              */
    	/*      n : 次数                                             */
    	/*      ct : 最大繰り返し回数                                */
    	/*      eps : 収束判定条件                                   */
    	/*      A : 対象とする行列                                   */
    	/*      A1, A2 : 作業域(nxnの行列),A1の対角要素が固有値   */
    	/*      X1, X2 : 作業域(nxnの行列),X1の各列が固有ベクトル */
    	/*      return : =0 : 正常                                   */
    	/*               =1 : 収束せず                               */
    	/*      coded by Y.Suganuma                                  */
    	/*************************************************************/
    	int Jacobi(int n, int ct, double eps, double A[][], double A1[][], double A2[][],
    	                  double X1[][], double X2[][])
    	{
    		double max, s, t, v, sn, cs;
    		int i1, i2, k = 0, ind = 1, p = 0, q = 0;
    					// 初期設定
    		for (i1 = 0; i1 < n; i1++) {
    			for (i2 = 0; i2 < n; i2++) {
    				A1[i1][i2] = A[i1][i2];
    				X1[i1][i2] = 0.0;
    			}
    			X1[i1][i1] = 1.0;
    		}
    					// 計算
    		while (ind > 0 && k < ct) {
    						// 最大要素の探索
    			max = 0.0;
    			for (i1 = 0; i1 < n; i1++) {
    				for (i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						if (Math.abs(A1[i1][i2]) > max) {
    							max = Math.abs(A1[i1][i2]);
    							p   = i1;
    							q   = i2;
    						}
    					}
    				}
    			}
    						// 収束判定
    							// 収束した
    			if (max < eps)
    				ind = 0;
    							// 収束しない
    			else {
    								// 準備
    				s  = -A1[p][q];
    				t  = 0.5 * (A1[p][p] - A1[q][q]);
    				v  = Math.abs(t) / Math.sqrt(s * s + t * t);
    				sn = Math.sqrt(0.5 * (1.0 - v));
    				if (s*t < 0.0)
    					sn = -sn;
    				cs = Math.sqrt(1.0 - sn * sn);
    								// Akの計算
    				for (i1 = 0; i1 < n; i1++) {
    					if (i1 == p) {
    						for (i2 = 0; i2 < n; i2++) {
    							if (i2 == p)
    								A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                               2.0 * A1[p][q] * sn * cs;
    							else if (i2 == q)
    								A2[p][q] = 0.0;
    							else
    								A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
    						}
    					}
    					else if (i1 == q) {
    						for (i2 = 0; i2 < n; i2++) {
    							if (i2 == q)
    								A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                               2.0 * A1[p][q] * sn * cs;
    							else if (i2 == p)
    								A2[q][p] = 0.0;
    							else
    								A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
    						}
    					}
    					else {
    						for (i2 = 0; i2 < n; i2++) {
    							if (i2 == p)
    								A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
    							else if (i2 == q)
    								A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
    							else
    								A2[i1][i2] = A1[i1][i2];
    						}
    					}
    				}
    								// Xkの計算
    				for (i1 = 0; i1 < n; i1++) {
    					for (i2 = 0; i2 < n; i2++) {
    						if (i2 == p)
    							X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
    						else if (i2 == q)
    							X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
    						else
    							X2[i1][i2] = X1[i1][i2];
    					}
    				}
    								// 次のステップへ
    				k++;
    				for (i1 = 0; i1 < n; i1++) {
    					for (i2 = 0; i2 < n; i2++) {
    						A1[i1][i2] = A2[i1][i2];
    						X1[i1][i2] = X2[i1][i2];
    					}
    				}
    			}
    		}
    
    		return ind;
    	}
    }
    			

  3. JavaScript

      ここをクリックすると,任意のデータに対して画面上で実行することができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>主成分分析</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		function main()
    		{
    					// データ
    			let p = parseInt(document.getElementById("p").value);
    			let n = parseInt(document.getElementById("n").value);
    			let r = new Array();
    			let x = new Array();
    			let a = new Array();
    			for (let i1 = 0; i1 < p; i1++) {
    				x[i1] = new Array();
    				a[i1] = new Array();
    			}
    			let ss = (document.getElementById("data").value).split(/ {1,}|\n{1,}/);
    			let k = 0;
    			for (let i1 = 0; i1 < n; i1++) {
    				for (let i2 = 0; i2 < p; i2++) {
    					x[i2][i1] = parseFloat(ss[k]);
    					k++;
    				}
    			}
    					// 計算
    			let sw = principal(p, n, x, r, a, 1.0e-10, 200);
    
    			if (sw > 0)
    				document.getElementById("ans").value = "収束しません";
    			else {
    				let str = "";
    				for (let i1 = 0; i1 < p; i1++) {
    					str = str + "主成分 " + r[i1] + "\n";
    					str = str + "  係数";
    					for (let i2 = 0; i2 < p; i2++)
    						str = str + " " + a[i1][i2];
    					str = str + "\n";
    				}
    				document.getElementById("ans").value = str;
    			}
    		}
    
    		/***********************************/
    		/* 主成分分析                      */
    		/*      p : 変数の数               */
    		/*      n : データの数             */
    		/*      x : データ                 */
    		/*      r : 分散(主成分)         */
    		/*      a : 係数                   */
    		/*      eps : 収束性を判定する規準 */
    		/*      ct : 最大繰り返し回数      */
    		/*      return : =0 : 正常         */
    		/*               =1 : エラー       */
    		/***********************************/
    		function principal(p, n, x, r, a, eps, ct)
    		{
    			let mean;
    			let s2;
    			let i1;
    			let i2;
    			let i3;
    			let sw = 0;
    						// 領域の確保
    			let C  = new Array();
    			let A1 = new Array();
    			let A2 = new Array();
    			let X1 = new Array();
    			let X2 = new Array();
    			for (let i1 = 0; i1 < p; i1++) {
    				C[i1]  = new Array();
    				A1[i1] = new Array();
    				A2[i1] = new Array();
    				X1[i1] = new Array();
    				X2[i1] = new Array();
    			}
    						// データの基準化
    			for (i1 = 0; i1 < p; i1++) {
    				mean = 0.0;
    				s2   = 0.0;
    				for (i2 = 0; i2 < n; i2++) {
    					mean += x[i1][i2];
    					s2   += x[i1][i2] * x[i1][i2];
    				}
    				mean /= n;
    				s2   /= n;
    				s2    = n * (s2 - mean * mean) / (n - 1);
    				s2    = Math.sqrt(s2);
    				for (i2 = 0; i2 < n; i2++)
    					x[i1][i2] = (x[i1][i2] - mean) / s2;
    			}
    						// 分散強分散行列の計算
    			for (i1 = 0; i1 < p; i1++) {
    				for (i2 = i1; i2 < p; i2++) {
    					s2 = 0.0;
    					for (i3 = 0; i3 < n; i3++)
    						s2 += x[i1][i3] * x[i2][i3];
    					s2        /= (n - 1);
    					C[i1][i2]  = s2;
    					if (i1 != i2)
    						C[i2][i1] = s2;
    				}
    			}
    						// 固有値と固有ベクトルの計算(ヤコビ法)
    			sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
    
    			if (sw == 0) {
    				for (i1 = 0; i1 < p; i1++) {
    					r[i1] = A1[i1][i1];
    					for (i2 = 0; i2 < p; i2++)
    						a[i1][i2] = X1[i2][i1];
    				}
    			}
    
    			return sw;
    		}
    
    		/*************************************************************/
    		/* 実対称行列の固有値・固有ベクトル(ヤコビ法)              */
    		/*      n : 次数                                             */
    		/*      ct : 最大繰り返し回数                                */
    		/*      eps : 収束判定条件                                   */
    		/*      A : 対象とする行列                                   */
    		/*      A1, A2 : 作業域(nxnの行列),A1の対角要素が固有値   */
    		/*      X1, X2 : 作業域(nxnの行列),X1の各列が固有ベクトル */
    		/*      return : =0 : 正常                                   */
    		/*               =1 : 収束せず                               */
    		/*      coded by Y.Suganuma                                  */
    		/*************************************************************/
    		function Jacobi(n, ct, eps, A, A1, A2, X1, X2)
    		{
    			let max;
    			let s;
    			let t;
    			let v;
    			let sn;
    			let cs;
    			let i1;
    			let i2;
    			let k = 0;
    			let ind = 1;
    			let p = 0;
    			let q = 0;
    						// 初期設定
    			for (i1 = 0; i1 < n; i1++) {
    				for (i2 = 0; i2 < n; i2++) {
    					A1[i1][i2] = A[i1][i2];
    					X1[i1][i2] = 0.0;
    				}
    				X1[i1][i1] = 1.0;
    			}
    						// 計算
    			while (ind > 0 && k < ct) {
    							// 最大要素の探索
    				max = 0.0;
    				for (i1 = 0; i1 < n; i1++) {
    					for (i2 = 0; i2 < n; i2++) {
    						if (i2 != i1) {
    							if (Math.abs(A1[i1][i2]) > max) {
    								max = Math.abs(A1[i1][i2]);
    								p   = i1;
    								q   = i2;
    							}
    						}
    					}
    				}
    							// 収束判定
    								// 収束した
    				if (max < eps)
    					ind = 0;
    								// 収束しない
    				else {
    									// 準備
    					s  = -A1[p][q];
    					t  = 0.5 * (A1[p][p] - A1[q][q]);
    					v  = Math.abs(t) / Math.sqrt(s * s + t * t);
    					sn = Math.sqrt(0.5 * (1.0 - v));
    					if (s*t < 0.0)
    						sn = -sn;
    					cs = Math.sqrt(1.0 - sn * sn);
    									// Akの計算
    					for (i1 = 0; i1 < n; i1++) {
    						if (i1 == p) {
    							for (i2 = 0; i2 < n; i2++) {
    								if (i2 == p)
    									A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                               	2.0 * A1[p][q] * sn * cs;
    								else if (i2 == q)
    									A2[p][q] = 0.0;
    								else
    									A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
    							}
    						}
    						else if (i1 == q) {
    							for (i2 = 0; i2 < n; i2++) {
    								if (i2 == q)
    									A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                               	2.0 * A1[p][q] * sn * cs;
    								else if (i2 == p)
    									A2[q][p] = 0.0;
    								else
    									A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
    							}
    						}
    						else {
    							for (i2 = 0; i2 < n; i2++) {
    								if (i2 == p)
    									A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
    								else if (i2 == q)
    									A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
    								else
    									A2[i1][i2] = A1[i1][i2];
    							}
    						}
    					}
    									// Xkの計算
    					for (i1 = 0; i1 < n; i1++) {
    						for (i2 = 0; i2 < n; i2++) {
    							if (i2 == p)
    								X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
    							else if (i2 == q)
    								X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
    							else
    								X2[i1][i2] = X1[i1][i2];
    						}
    					}
    									// 次のステップへ
    					k++;
    					for (i1 = 0; i1 < n; i1++) {
    						for (i2 = 0; i2 < n; i2++) {
    							A1[i1][i2] = A2[i1][i2];
    							X1[i1][i2] = X2[i1][i2];
    						}
    					}
    				}
    			}
    
    			return ind;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>主成分分析</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,テキストエリアに与えられた 100 個のデータに対して主成分分析を行う場合に対する値が設定されています.他の問題を実行する場合は,それらを適切に修正してください.
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		変数の数:<INPUT ID="p" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="4"> 
    		データの数:<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="100"><BR><BR>
    		データ:<TEXTAREA ID="data" COLS="30" ROWS="15" STYLE="font-size: 100%">
    66 22 44 31
    25 74 17 81
    50 23 53 71
    25 57 19 81
    74 47 64 47
    39 33 48 46
    14 22 9 69
    67 60 49 26
    42 40 77 65
    11 80 0 86
    32 0 43 74
    68 69 44 68
    24 49 9 71
    42 74 28 46
    60 58 73 28
    36 37 33 68
    24 44 19 83
    30 40 31 50
    55 40 60 49
    63 47 94 41
    72 30 100 45
    19 22 13 75
    43 39 43 34
    90 83 92 31
    51 77 52 82
    53 70 34 31
    28 51 53 44
    40 62 42 79
    31 48 22 68
    57 29 51 30
    64 89 57 42
    49 82 72 29
    53 31 55 43
    79 52 70 10
    45 19 43 57
    35 34 34 89
    4 69 0 100
    49 49 66 66
    92 82 97 6
    5 89 0 100
    65 26 83 28
    56 36 64 38
    48 50 25 22
    30 30 15 55
    40 65 38 42
    14 67 9 67
    84 96 90 8
    53 64 51 54
    50 89 60 52
    76 41 68 9
    49 40 53 49
    78 66 66 17
    76 58 90 29
    41 15 40 49
    63 60 55 33
    40 36 49 67
    78 54 71 18
    62 72 69 12
    64 47 42 53
    56 64 9 15
    77 35 56 25
    44 12 46 87
    80 9 56 19
    36 21 52 78
    48 63 64 48
    43 61 50 47
    58 23 28 50
    90 12 100 0
    13 33 11 77
    67 44 48 28
    75 45 68 17
    81 22 89 9
    46 45 59 55
    56 49 64 55
    65 62 72 27
    34 49 29 77
    45 33 60 63
    20 45 14 99
    33 38 26 87
    44 51 69 52
    64 57 64 48
    44 64 51 28
    63 48 56 11
    29 39 33 84
    40 48 51 54
    40 38 26 62
    68 46 61 26
    58 45 68 48
    64 44 77 63
    59 62 44 66
    81 53 93 19
    23 34 12 68
    51 35 55 46
    74 70 84 17
    42 33 56 44
    46 31 46 53
    33 57 38 63
    40 24 20 42
    53 36 60 31
    0 34 0 100
    		</TEXTAREA> 
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
    		<TEXTAREA ID="ans" COLS="50" ROWS="5" STYLE="font-size: 100%;"></TEXTAREA>
    	</DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /****************************/
    /* 主成分分析               */
    /*      coded by Y.Suganuma */
    /****************************/
    
    	fscanf(STDIN, "%d %d", $p, $n);   // 変数の数とデータの数
    
    	$r = array($p);
    	$x = array($p);
    	$a = array($p);
    	for ($i1 = 0; $i1 < $p; $i1++) {
    		$x[$i1]  = array($n);
    		$a[$i1]  = array($p);
    	}
    
    	for ($i1 = 0; $i1 < $n; $i1++) {   // データ
    		$str       = trim(fgets(STDIN));
    		$x[0][$i1] = floatval(strtok($str, " "));
    		for ($i2 = 1; $i2 < $p; $i2++)
    			$x[$i2][$i1] = floatval(strtok(" "));
    	}
    
    	$sw = principal($p, $n, $x, $r, $a, 1.0e-10, 200);
    
    	if ($sw == 0) {
    		for ($i1 = 0; $i1 < $p; $i1++) {
    			printf("主成分 %f", $r[$i1]);
    			printf(" 係数");
    			for ($i2 = 0; $i2 < $p; $i2++)
    				printf(" %f", $a[$i1][$i2]);
    			printf("\n");
    		}
    	}
    	else
    		printf("***error  解を求めることができませんでした\n");
    
    /***********************************/
    /* 主成分分析                      */
    /*      p : 変数の数               */
    /*      n : データの数             */
    /*      x : データ                 */
    /*      r : 分散(主成分)         */
    /*      a : 係数                   */
    /*      eps : 収束性を判定する規準 */
    /*      ct : 最大繰り返し回数      */
    /*      return : =0 : 正常         */
    /*               =1 : エラー       */
    /***********************************/
    function principal($p, $n,  $x, &$r, &$a, $eps, $ct)
    {
    	$sw = 0;
    					// 領域の確保
    	$C  = array($p);
    	$A1 = array($p);
    	$A2 = array($p);
    	$X1 = array($p);
    	$X2 = array($p);
    	for ($i1 = 0; $i1 < $p; $i1++) {
    		$C[$i1]  = array($p);
    		$A1[$i1] = array($p);
    		$A2[$i1] = array($p);
    		$X1[$i1] = array($p);
    		$X2[$i1] = array($p);
    	}
    					// データの基準化
    	for ($i1 = 0; $i1 < $p; $i1++) {
    		$mean = 0.0;
    		$s2   = 0.0;
    		for ($i2 = 0; $i2 < $n; $i2++) {
    			$mean += $x[$i1][$i2];
    			$s2   += $x[$i1][$i2] * $x[$i1][$i2];
    		}
    		$mean /= $n;
    		$s2   /= $n;
    		$s2    = $n * ($s2 - $mean * $mean) / ($n - 1);
    		$s2    = sqrt($s2);
    		for ($i2 = 0; $i2 < $n; $i2++)
    			$x[$i1][$i2] = ($x[$i1][$i2] - $mean) / $s2;
    	}
    					// 分散強分散行列の計算
    	for ($i1 = 0; $i1 < $p; $i1++) {
    		for ($i2 = $i1; $i2 < $p; $i2++) {
    			$s2 = 0.0;
    			for ($i3 = 0; $i3 < $n; $i3++)
    				$s2 += $x[$i1][$i3] * $x[$i2][$i3];
    			$s2          /= ($n - 1);
    			$C[$i1][$i2]  = $s2;
    			if ($i1 != $i2)
    				$C[$i2][$i1] = $s2;
    		}
    	}
    					// 固有値と固有ベクトルの計算(ヤコビ法)
    	$sw = Jacobi($p, $ct, $eps, $C, $A1, $A2, $X1, $X2);
    
    	if ($sw == 0) {
    		for ($i1 = 0; $i1 < $p; $i1++) {
    			$r[$i1] = $A1[$i1][$i1];
    			for ($i2 = 0; $i2 < $p; $i2++)
    				$a[$i1][$i2] = $X1[$i2][$i1];
    		}
    	}
    
    	return $sw;
    }
    
    /*************************************************************/
    /* 実対称行列の固有値・固有ベクトル(ヤコビ法)              */
    /*      n : 次数                                             */
    /*      ct : 最大繰り返し回数                                */
    /*      eps : 収束判定条件                                   */
    /*      A : 対象とする行列                                   */
    /*      A1, A2 : 作業域(nxnの行列),A1の対角要素が固有値   */
    /*      X1, X2 : 作業域(nxnの行列),X1の各列が固有ベクトル */
    /*      return : =0 : 正常                                   */
    /*               =1 : 収束せず                               */
    /*      coded by Y.Suganuma                                  */
    /*************************************************************/
    function Jacobi($n, $ct, $eps, $A, &$A1, $A2, &$X1, $X2)
    {
    					// 初期設定
    	$k   = 0;
    	$ind = 1;
    	$p   = 0;
    	$q   = 0;
    
    	for ($i1 = 0; $i1 < $n; $i1++) {
    		for ($i2 = 0; $i2 < $n; $i2++) {
    			$A1[$i1][$i2] = $A[$i1][$i2];
    			$X1[$i1][$i2] = 0.0;
    		}
    		$X1[$i1][$i1] = 1.0; 
    	}
    					// 計算
    	while ($ind > 0 && $k < $ct) {
    						// 最大要素の探索
    		$max = 0.0;
    		for ($i1 = 0; $i1 < $n; $i1++) {
    			for ($i2 = 0; $i2 < $n; $i2++) {
    				if ($i2 != $i1) {
    					if (abs($A1[$i1][$i2]) > $max) {
    						$max = abs($A1[$i1][$i2]);
    						$p   = $i1;
    						$q   = $i2;
    					}
    				}
    			}
    		}
    						// 収束判定
    							// 収束した
    		if ($max < $eps)
    			$ind = 0;
    							// 収束しない
    		else {
    								// 準備
    			$s  = -$A1[$p][$q];
    			$t  = 0.5 * ($A1[$p][$p] - $A1[$q][$q]);
    			$v  = abs($t) / sqrt($s * $s + $t * $t);
    			$sn = sqrt(0.5 * (1.0 - $v));
    			if ($s*$t < 0.0)
    				$sn = -$sn;
    			$cs = sqrt(1.0 - $sn * $sn);
    								// Akの計算
    			for ($i1 = 0; $i1 < $n; $i1++) {
    				if ($i1 == $p) {
    					for ($i2 = 0; $i2 < $n; $i2++) {
    						if ($i2 == $p)
    							$A2[$p][$p] = $A1[$p][$p] * $cs * $cs + $A1[$q][$q] * $sn * $sn -
                                              2.0 * $A1[$p][$q] * $sn * $cs;
    						else if ($i2 == $q)
    							$A2[$p][$q] = 0.0;
    						else
    							$A2[$p][$i2] = $A1[$p][$i2] * $cs - $A1[$q][$i2] * $sn;
    					}
    				}
    				else if ($i1 == $q) {
    					for ($i2 = 0; $i2 < $n; $i2++) {
    						if ($i2 == $q)
    							$A2[$q][$q] = $A1[$p][$p] * $sn * $sn + $A1[$q][$q] * $cs * $cs +
                                              2.0 * $A1[$p][$q] * $sn * $cs;
    						else if ($i2 == $p)
    							$A2[$q][$p] = 0.0;
    						else
    							$A2[$q][$i2] = $A1[$q][$i2] * $cs + $A1[$p][$i2] * $sn;
    					}
    				}
    				else {
    					for ($i2 = 0; $i2 < $n; $i2++) {
    						if ($i2 == $p)
    							$A2[$i1][$p] = $A1[$i1][$p] * $cs - $A1[$i1][$q] * $sn;
    						else if ($i2 == $q)
    							$A2[$i1][$q] = $A1[$i1][$q] * $cs + $A1[$i1][$p] * $sn;
    						else
    							$A2[$i1][$i2] = $A1[$i1][$i2];
    					}
    				}
    			}
    								// Xkの計算
    			for ($i1 = 0; $i1 < $n; $i1++) {
    				for ($i2 = 0; $i2 < $n; $i2++) {
    					if ($i2 == $p)
    						$X2[$i1][$p] = $X1[$i1][$p] * $cs - $X1[$i1][$q] * $sn;
    					else if ($i2 == $q)
    						$X2[$i1][$q] = $X1[$i1][$q] * $cs + $X1[$i1][$p] * $sn;
    					else
    						$X2[$i1][$i2] = $X1[$i1][$i2];
    				}
    			}
    								// 次のステップへ
    			$k++;
    			for ($i1 = 0; $i1 < $n; $i1++) {
    				for ($i2 = 0; $i2 < $n; $i2++) {
    					$A1[$i1][$i2] = $A2[$i1][$i2];
    					$X1[$i1][$i2] = $X2[$i1][$i2];
    				}
    			}
    		}
    	}
    
    	return $ind;
    }
    
    /*
    ---------データ例(コメント部分を除いて下さい)---------
    4 100   // 変数の数(p)とデータの数(n)
    66 22 44 31   // x1, x2, x3, x4
    25 74 17 81
    50 23 53 71
    25 57 19 81
    74 47 64 47
    39 33 48 46
    14 22 9 69
    67 60 49 26
    42 40 77 65
    11 80 0 86
    32 0 43 74
    68 69 44 68
    24 49 9 71
    42 74 28 46
    60 58 73 28
    36 37 33 68
    24 44 19 83
    30 40 31 50
    55 40 60 49
    63 47 94 41
    72 30 100 45
    19 22 13 75
    43 39 43 34
    90 83 92 31
    51 77 52 82
    53 70 34 31
    28 51 53 44
    40 62 42 79
    31 48 22 68
    57 29 51 30
    64 89 57 42
    49 82 72 29
    53 31 55 43
    79 52 70 10
    45 19 43 57
    35 34 34 89
    4 69 0 100
    49 49 66 66
    92 82 97 6
    5 89 0 100
    65 26 83 28
    56 36 64 38
    48 50 25 22
    30 30 15 55
    40 65 38 42
    14 67 9 67
    84 96 90 8
    53 64 51 54
    50 89 60 52
    76 41 68 9
    49 40 53 49
    78 66 66 17
    76 58 90 29
    41 15 40 49
    63 60 55 33
    40 36 49 67
    78 54 71 18
    62 72 69 12
    64 47 42 53
    56 64 9 15
    77 35 56 25
    44 12 46 87
    80 9 56 19
    36 21 52 78
    48 63 64 48
    43 61 50 47
    58 23 28 50
    90 12 100 0
    13 33 11 77
    67 44 48 28
    75 45 68 17
    81 22 89 9
    46 45 59 55
    56 49 64 55
    65 62 72 27
    34 49 29 77
    45 33 60 63
    20 45 14 99
    33 38 26 87
    44 51 69 52
    64 57 64 48
    44 64 51 28
    63 48 56 11
    29 39 33 84
    40 48 51 54
    40 38 26 62
    68 46 61 26
    58 45 68 48
    64 44 77 63
    59 62 44 66
    81 53 93 19
    23 34 12 68
    51 35 55 46
    74 70 84 17
    42 33 56 44
    46 31 46 53
    33 57 38 63
    40 24 20 42
    53 36 60 31
    0 34 0 100
    */
    
    ?>
    			

  5. Ruby

    ############################
    # 主成分分析
    #      coded by Y.Suganuma
    ############################
    
    #************************************************************/
    # 実対称行列の固有値・固有ベクトル(ヤコビ法)              */
    #      n : 次数                                             */
    #      ct : 最大繰り返し回数                                */
    #      eps : 収束判定条件                                   */
    #      a : 対象とする行列                                   */
    #      a1, a2 : 作業域(nxnの行列),a1の対角要素が固有値   */
    #      x1, x2 : 作業域(nxnの行列),x1の各列が固有ベクトル */
    #      return : =0 : 正常                                   */
    #               =1 : 収束せず                               */
    #      coded by Y.Suganuma                                  */
    #************************************************************/
    def Jacobi(n, ct, eps, a, a1, a2, x1, x2)
    					# 初期設定
    	k   = 0
    	ind = 1
    	p   = 0
    	q   = 0
    	for i1 in 0 ... n
    		for i2 in 0 ... n
    			a1[i1][i2] = a[i1][i2]
    			x1[i1][i2] = 0.0
    		end
    		x1[i1][i1] = 1.0
    	end
    					# 計算
    	while ind > 0 && k < ct
    						# 最大要素の探索
    		max = 0.0
    		for i1 in 0 ... n
    			for i2 in 0 ... n
    				if i2 != i1
    					if a1[i1][i2].abs() > max
    						max = a1[i1][i2].abs()
    						p   = i1
    						q   = i2
    					end
    				end
    			end
    		end
    						# 収束判定
    							# 収束した
    		if max < eps
    			ind = 0
    							# 収束しない
    		else
    								# 準備
    			s  = -a1[p][q]
    			t  = 0.5 * (a1[p][p] - a1[q][q])
    			v  = t.abs() / Math.sqrt(s * s + t * t)
    			sn = Math.sqrt(0.5 * (1.0 - v))
    			if s*t < 0.0
    				sn = -sn
    			end
    			cs = Math.sqrt(1.0 - sn * sn)
    								# akの計算
    			for i1 in 0 ... n
    				if i1 == p
    					for i2 in 0 ... n
    						if i2 == p
    							a2[p][p] = a1[p][p] * cs * cs + a1[q][q] * sn * sn -
                                           2.0 * a1[p][q] * sn * cs
    						elsif i2 == q
    							a2[p][q] = 0.0
    						else
    							a2[p][i2] = a1[p][i2] * cs - a1[q][i2] * sn
    						end
    					end
    				elsif i1 == q
    					for i2 in 0 ... n
    						if (i2 == q)
    							a2[q][q] = a1[p][p] * sn * sn + a1[q][q] * cs * cs +
                                           2.0 * a1[p][q] * sn * cs
    						elsif i2 == p
    							a2[q][p] = 0.0
    						else
    							a2[q][i2] = a1[q][i2] * cs + a1[p][i2] * sn
    						end
    					end
    				else
    					for i2 in 0 ... n
    						if i2 == p
    							a2[i1][p] = a1[i1][p] * cs - a1[i1][q] * sn
    						elsif i2 == q
    							a2[i1][q] = a1[i1][q] * cs + a1[i1][p] * sn
    						else
    							a2[i1][i2] = a1[i1][i2]
    						end
    					end
    				end
    			end
    								# xkの計算
    			for i1 in 0 ... n
    				for i2 in 0 ... n
    					if i2 == p
    						x2[i1][p] = x1[i1][p] * cs - x1[i1][q] * sn
    					elsif i2 == q
    						x2[i1][q] = x1[i1][q] * cs + x1[i1][p] * sn
    					else
    						x2[i1][i2] = x1[i1][i2]
    					end
    				end
    			end
    								# 次のステップへ
    			k += 1
    			for i1 in 0 ... n
    				for i2 in 0 ... n
    					a1[i1][i2] = a2[i1][i2]
    					x1[i1][i2] = x2[i1][i2]
    				end
    			end
    		end
    	end
    
    	return ind
    end
    
    ###################################
    # 主成分分析
    #      p : 変数の数
    #      n : データの数
    #      x : データ
    #      r : 分散(主成分)
    #      a : 係数
    #      eps : 収束性を判定する規準
    #      ct : 最大繰り返し回数
    #      return : =0 : 正常
    #               =1 : エラー
    #      coded by Y.Suganuma
    ###################################
    def principal(p, n, x, r, a, eps, ct)
    
    	sw = 0
    			# 領域の確保
    	c  = Array.new(p)
    	a1 = Array.new(p)
    	a2 = Array.new(p)
    	x1 = Array.new(p)
    	x2 = Array.new(p)
    	for i1 in 0 ... p
    		c[i1]  = Array.new(p)
    		a1[i1] = Array.new(p)
    		a2[i1] = Array.new(p)
    		x1[i1] = Array.new(p)
    		x2[i1] = Array.new(p)
    	end
    			# データの基準化
    	for i1 in 0 ... p
    		mean = 0.0
    		s2   = 0.0
    		for i2 in 0 ... n
    			mean += x[i1][i2]
    			s2   += x[i1][i2] * x[i1][i2]
    		end
    		mean /= n
    		s2   /= n
    		s2    = n * (s2 - mean * mean) / (n - 1)
    		s2    = Math.sqrt(s2)
    		for i2 in 0 ... n
    			x[i1][i2] = (x[i1][i2] - mean) / s2
    		end
    	end
    			# 分散強分散行列の計算
    	for i1 in 0 ... p
    		for i2 in 0 ... p
    			s2 = 0.0
    			for i3 in 0 ... n
    				s2 += x[i1][i3] * x[i2][i3]
    			end
    			s2        /= (n - 1)
    			c[i1][i2]  = s2
    			if i1 != i2
    				c[i2][i1] = s2
    			end
    		end
    	end
    			# 固有値と固有ベクトルの計算(ヤコビ法)
    	sw = Jacobi(p, ct, eps, c, a1, a2, x1, x2)
    
    	if sw == 0
    		for i1 in 0 ... p
    			r[i1] = a1[i1][i1]
    			for i2 in 0 ... p
    				a[i1][i2] = x1[i2][i1]
    			end
    		end
    	end
    
    	return sw
    end
    
    ss = gets().split(" ")
    p  = Integer(ss[0])   # 変数の数
    n  = Integer(ss[1])   # データの数
    
    r  = Array.new(p)
    x  = Array.new(p)
    a  = Array.new(p)
    for i1 in 0 ... p
    	x[i1] = Array.new(n)
    	a[i1] = Array.new(p)
    end
    
    for i1 in 0 ... n   # データ
    	ss = gets().split(" ")
    	for i2 in 0 ... p
    		x[i2][i1] = Float(ss[i2])
    	end
    end
    
    sw = principal(p, n, x, r, a, 1.0e-10, 200)
    
    if sw == 0
    	for i1 in 0 ... p
    		print("主成分 " + String(r[i1]))
    		print(" 係数")
    		for i2 in 0 ... p
    			print(" " + String(a[i1][i2]))
    		end
    		print("\n")
    	end
    else
    	print("***error  解を求めることができませんでした")
    end
    
    =begin
    ---------データ例(コメント部分を除いて下さい)---------
    4 100   # 変数の数(p)とデータの数(n)
    66 22 44 31   # x1, x2, x3, x4
    25 74 17 81
    50 23 53 71
    25 57 19 81
    74 47 64 47
    39 33 48 46
    14 22 9 69
    67 60 49 26
    42 40 77 65
    11 80 0 86
    32 0 43 74
    68 69 44 68
    24 49 9 71
    42 74 28 46
    60 58 73 28
    36 37 33 68
    24 44 19 83
    30 40 31 50
    55 40 60 49
    63 47 94 41
    72 30 100 45
    19 22 13 75
    43 39 43 34
    90 83 92 31
    51 77 52 82
    53 70 34 31
    28 51 53 44
    40 62 42 79
    31 48 22 68
    57 29 51 30
    64 89 57 42
    49 82 72 29
    53 31 55 43
    79 52 70 10
    45 19 43 57
    35 34 34 89
    4 69 0 100
    49 49 66 66
    92 82 97 6
    5 89 0 100
    65 26 83 28
    56 36 64 38
    48 50 25 22
    30 30 15 55
    40 65 38 42
    14 67 9 67
    84 96 90 8
    53 64 51 54
    50 89 60 52
    76 41 68 9
    49 40 53 49
    78 66 66 17
    76 58 90 29
    41 15 40 49
    63 60 55 33
    40 36 49 67
    78 54 71 18
    62 72 69 12
    64 47 42 53
    56 64 9 15
    77 35 56 25
    44 12 46 87
    80 9 56 19
    36 21 52 78
    48 63 64 48
    43 61 50 47
    58 23 28 50
    90 12 100 0
    13 33 11 77
    67 44 48 28
    75 45 68 17
    81 22 89 9
    46 45 59 55
    56 49 64 55
    65 62 72 27
    34 49 29 77
    45 33 60 63
    20 45 14 99
    33 38 26 87
    44 51 69 52
    64 57 64 48
    44 64 51 28
    63 48 56 11
    29 39 33 84
    40 48 51 54
    40 38 26 62
    68 46 61 26
    58 45 68 48
    64 44 77 63
    59 62 44 66
    81 53 93 19
    23 34 12 68
    51 35 55 46
    74 70 84 17
    42 33 56 44
    46 31 46 53
    33 57 38 63
    40 24 20 42
    53 36 60 31
    0 34 0 100
    =end
    			

  6. Python

    # -*- coding: UTF-8 -*-
    import numpy as np
    import sys
    from math import *
    
    ############################################
    # 実対称行列の固有値・固有ベクトル(ヤコビ法)
    #      n : 次数
    #      ct : 最大繰り返し回数
    #      eps : 収束判定条件
    #      A : 対象とする行列
    #      A1, A2 : 作業域(nxnの行列),A1の対角要素が固有値
    #      X1, X2 : 作業域(nxnの行列),X1の各列が固有ベクトル
    #      return : =0 : 正常
    #               =1 : 収束せず
    #      coded by Y.Suganuma
    ############################################
    
    def Jacobi(n, ct, eps, A, A1, A2, X1, X2) :
    			# 初期設定
    	k   = 0
    	ind = 1
    	p   = 0
    	q   = 0
    	for i1 in range(0, n) :
    		for i2 in range(0, n) :
    			A1[i1][i2] = A[i1][i2]
    			X1[i1][i2] = 0.0
    		X1[i1][i1] = 1.0
    			# 計算
    	while ind > 0 and k < ct :
    				# 最大要素の探索
    		max = 0.0
    		for i1 in range(0, n) :
    			for i2 in range(0, n) :
    				if i2 != i1 :
    					if abs(A1[i1][i2]) > max :
    						max = fabs(A1[i1][i2])
    						p   = i1
    						q   = i2
    				# 収束判定
    					# 収束した
    		if max < eps :
    			ind = 0
    					# 収束しない
    		else :
    						# 準備
    			s  = -A1[p][q]
    			t  = 0.5 * (A1[p][p] - A1[q][q])
    			v  = fabs(t) / sqrt(s * s + t * t)
    			sn = sqrt(0.5 * (1.0 - v))
    			if s*t < 0.0 :
    				sn = -sn
    			cs = sqrt(1.0 - sn * sn)
    						# Akの計算
    			for i1 in range(0, n) :
    				if i1 == p :
    					for i2 in range(0, n) :
    						if i2 == p :
    							A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn - 2.0 * A1[p][q] * sn * cs
    						elif i2 == q :
    							A2[p][q] = 0.0
    						else :
    							A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn
    				elif i1 == q :
    					for i2 in range(0, n) :
    						if i2 == q :
    							A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs + 2.0 * A1[p][q] * sn * cs
    						elif i2 == p :
    							A2[q][p] = 0.0
    						else :
    							A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn
    				else :
    					for i2 in range(0, n) :
    						if i2 == p :
    							A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn
    						elif i2 == q :
    							A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn
    						else :
    							A2[i1][i2] = A1[i1][i2]
    						# Xkの計算
    			for i1 in range(0, n) :
    				for i2 in range(0, n) :
    					if i2 == p :
    						X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn
    					elif i2 == q :
    						X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn
    					else :
    						X2[i1][i2] = X1[i1][i2]
    						# 次のステップへ
    			k += 1
    			for i1 in range(0, n) :
    				for i2 in range(0, n) :
    					A1[i1][i2] = A2[i1][i2]
    					X1[i1][i2] = X2[i1][i2]
    
    	return ind
    
    ###################################
    # 主成分分析
    #      p : 変数の数
    #      n : データの数
    #      x : データ
    #      r : 分散(主成分)
    #      a : 係数
    #      eps : 収束性を判定する規準
    #      ct : 最大繰り返し回数
    #      return : =0 : 正常
    #               =1 : エラー
    #      coded by Y.Suganuma
    ###################################
    def principal(p, n, x, r, a, eps, ct) :
    
    	sw = 0
    			# 領域の確保
    	C  = np.empty((p, p), np.float)
    	A1 = np.empty((p, p), np.float)
    	A2 = np.empty((p, p), np.float)
    	X1 = np.empty((p, p), np.float)
    	X2 = np.empty((p, p), np.float)
    			# データの基準化
    	for i1 in range(0, p) :
    		mean = 0.0
    		s2   = 0.0
    		for i2 in range(0, n) :
    			mean += x[i1][i2]
    			s2   += x[i1][i2] * x[i1][i2]
    		mean /= n
    		s2   /= n
    		s2    = n * (s2 - mean * mean) / (n - 1)
    		s2    = sqrt(s2)
    		for i2 in range(0, n) :
    			x[i1][i2] = (x[i1][i2] - mean) / s2
    			# 分散強分散行列の計算
    	for i1 in range(0, p) :
    		for i2 in range(0, p) :
    			s2 = 0.0
    			for i3 in range(0, n) :
    				s2 += x[i1][i3] * x[i2][i3]
    			s2        /= (n - 1)
    			C[i1][i2]  = s2
    			if i1 != i2 :
    				C[i2][i1] = s2
    			# 固有値と固有ベクトルの計算(ヤコビ法)
    	sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2)
    
    	if sw == 0 :
    		for i1 in range(0, p) :
    			r[i1] = A1[i1][i1]
    			for i2 in range(0, p) :
    				a[i1][i2] = X1[i2][i1]
    
    	return sw
    
    ############################
    # 主成分分析
    #      coded by Y.Suganuma
    ############################
    
    line = sys.stdin.readline()
    ss   = line.split()
    p    = int(ss[0])   # 変数の数
    n    = int(ss[1])   # データの数
    
    r    = np.empty(p, np.float)
    x    = np.empty((p, n), np.float)
    a    = np.empty((p, p), np.float)
    
    for i1 in range(0, n) :   # データ
    	line = sys.stdin.readline()
    	ss   = line.split()
    	for i2 in range(0, p) :
    		x[i2][i1] = float(ss[i2])
    
    sw = principal(p, n, x, r, a, 1.0e-10, 200)
    
    if sw == 0 :
    	for i1 in range(0, p) :
    		print("主成分 " + str(r[i1]), end="")
    		print(" 係数", end="")
    		for i2 in range(0, p) :
    			print(" " + str(a[i1][i2]), end="")
    		print()
    else :
    	print("***error  解を求めることができませんでした")
    
    """
    ---------データ例(コメント部分を除いて下さい)---------
    4 100   # 変数の数(p)とデータの数(n)
    66 22 44 31   # x1, x2, x3, x4
    25 74 17 81
    50 23 53 71
    25 57 19 81
    74 47 64 47
    39 33 48 46
    14 22 9 69
    67 60 49 26
    42 40 77 65
    11 80 0 86
    32 0 43 74
    68 69 44 68
    24 49 9 71
    42 74 28 46
    60 58 73 28
    36 37 33 68
    24 44 19 83
    30 40 31 50
    55 40 60 49
    63 47 94 41
    72 30 100 45
    19 22 13 75
    43 39 43 34
    90 83 92 31
    51 77 52 82
    53 70 34 31
    28 51 53 44
    40 62 42 79
    31 48 22 68
    57 29 51 30
    64 89 57 42
    49 82 72 29
    53 31 55 43
    79 52 70 10
    45 19 43 57
    35 34 34 89
    4 69 0 100
    49 49 66 66
    92 82 97 6
    5 89 0 100
    65 26 83 28
    56 36 64 38
    48 50 25 22
    30 30 15 55
    40 65 38 42
    14 67 9 67
    84 96 90 8
    53 64 51 54
    50 89 60 52
    76 41 68 9
    49 40 53 49
    78 66 66 17
    76 58 90 29
    41 15 40 49
    63 60 55 33
    40 36 49 67
    78 54 71 18
    62 72 69 12
    64 47 42 53
    56 64 9 15
    77 35 56 25
    44 12 46 87
    80 9 56 19
    36 21 52 78
    48 63 64 48
    43 61 50 47
    58 23 28 50
    90 12 100 0
    13 33 11 77
    67 44 48 28
    75 45 68 17
    81 22 89 9
    46 45 59 55
    56 49 64 55
    65 62 72 27
    34 49 29 77
    45 33 60 63
    20 45 14 99
    33 38 26 87
    44 51 69 52
    64 57 64 48
    44 64 51 28
    63 48 56 11
    29 39 33 84
    40 48 51 54
    40 38 26 62
    68 46 61 26
    58 45 68 48
    64 44 77 63
    59 62 44 66
    81 53 93 19
    23 34 12 68
    51 35 55 46
    74 70 84 17
    42 33 56 44
    46 31 46 53
    33 57 38 63
    40 24 20 42
    53 36 60 31
    0 34 0 100
    """
    			

  7. C#

    /****************************/
    /* 主成分分析               */
    /*      coded by Y.Suganuma */
    /****************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	public Test1()
    	{
    		char[] charSep = new char[] {' '};
    					// 変数の数とデータの数
    		String[] str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
    		int p        = int.Parse(str[0]);
    		int n        = int.Parse(str[1]);
    
    		double[] r   = new double [p];
    		double[][] x = new double [p][];
    		double[][] a = new double [p][];
    		for (int i1 = 0; i1 < p; i1++) {
    			x[i1] = new double [n];
    			a[i1] = new double [p];
    		}
    					// データ
    		for (int i1 = 0; i1 < n; i1++) {
    			str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
    			for (int i2 = 0; i2 < p; i2++)
    				x[i2][i1] = double.Parse(str[i2]);
    		}
    
    		int sw = principal(p, n, x, r, a, 1.0e-10, 200);
    
    		if (sw == 0) {
    			for (int i1 = 0; i1 < p; i1++) {
    				Console.Write("主成分 " + r[i1]);
    				Console.Write(" 係数");
    				for (int i2 = 0; i2 < p; i2++)
    					Console.Write(" " + a[i1][i2]);
    				Console.WriteLine();
    			}
    		}
    		else
    			Console.WriteLine("***error  解を求めることができませんでした");
    	}
    
    	/***********************************/
    	/* 主成分分析                      */
    	/*      p : 変数の数               */
    	/*      n : データの数             */
    	/*      x : データ                 */
    	/*      r : 分散(主成分)         */
    	/*      a : 係数                   */
    	/*      eps : 収束性を判定する規準 */
    	/*      ct : 最大繰り返し回数      */
    	/*      return : =0 : 正常         */
    	/*               =1 : エラー       */
    	/***********************************/
    	int principal(int p, int n, double[][] x, double[] r, double[][] a, double eps, int ct)
    	{
    					// 領域の確保
    		double[][] C  = new double [p][];
    		double[][] A1 = new double [p][];
    		double[][] A2 = new double [p][];
    		double[][] X1 = new double [p][];
    		double[][] X2 = new double [p][];
    		for (int i1 = 0; i1 < p; i1++) {
    			C[i1]  = new double [p];
    			A1[i1] = new double [p];
    			A2[i1] = new double [p];
    			X1[i1] = new double [p];
    			X2[i1] = new double [p];
    		}
    					// データの基準化
    		for (int i1 = 0; i1 < p; i1++) {
    			double mean = 0.0;
    			double s2   = 0.0;
    			for (int i2 = 0; i2 < n; i2++) {
    				mean += x[i1][i2];
    				s2   += x[i1][i2] * x[i1][i2];
    			}
    			mean /= n;
    			s2   /= n;
    			s2    = n * (s2 - mean * mean) / (n - 1);
    			s2    = Math.Sqrt(s2);
    			for (int i2 = 0; i2 < n; i2++)
    				x[i1][i2] = (x[i1][i2] - mean) / s2;
    		}
    					// 分散強分散行列の計算
    		for (int i1 = 0; i1 < p; i1++) {
    			for (int i2 = i1; i2 < p; i2++) {
    				double s2 = 0.0;
    				for (int i3 = 0; i3 < n; i3++)
    					s2 += x[i1][i3] * x[i2][i3];
    				s2        /= (n - 1);
    				C[i1][i2]  = s2;
    				if (i1 != i2)
    					C[i2][i1] = s2;
    			}
    		}
    					// 固有値と固有ベクトルの計算(ヤコビ法)
    		int sw = Jacobi(p, ct, eps, C, A1, A2, X1, X2);
    
    		if (sw == 0) {
    			for (int i1 = 0; i1 < p; i1++) {
    				r[i1] = A1[i1][i1];
    				for (int i2 = 0; i2 < p; i2++)
    					a[i1][i2] = X1[i2][i1];
    			}
    		}
    
    		return sw;
    	}
    
    	/*************************************************************/
    	/* 実対称行列の固有値・固有ベクトル(ヤコビ法)              */
    	/*      n : 次数                                             */
    	/*      ct : 最大繰り返し回数                                */
    	/*      eps : 収束判定条件                                   */
    	/*      A : 対象とする行列                                   */
    	/*      A1, A2 : 作業域(nxnの行列),A1の対角要素が固有値   */
    	/*      X1, X2 : 作業域(nxnの行列),X1の各列が固有ベクトル */
    	/*      return : =0 : 正常                                   */
    	/*               =1 : 収束せず                               */
    	/*      coded by Y.Suganuma                                  */
    	/*************************************************************/
    	int Jacobi(int n, int ct, double eps, double[][] A, double[][] A1, double[][] A2,
    	           double[][] X1, double[][] X2)
    	{
    					// 初期設定
    		for (int i1 = 0; i1 < n; i1++) {
    			for (int i2 = 0; i2 < n; i2++) {
    				A1[i1][i2] = A[i1][i2];
    				X1[i1][i2] = 0.0;
    			}
    			X1[i1][i1] = 1.0;
    		}
    					// 計算
    		int k = 0, ind = 1;
    		while (ind > 0 && k < ct) {
    						// 最大要素の探索
    			double max = 0.0;
    			int p = 0, q = 0;
    			for (int i1 = 0; i1 < n; i1++) {
    				for (int i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						if (Math.Abs(A1[i1][i2]) > max) {
    							max = Math.Abs(A1[i1][i2]);
    							p   = i1;
    							q   = i2;
    						}
    					}
    				}
    			}
    						// 収束判定
    							// 収束した
    			if (max < eps)
    				ind = 0;
    							// 収束しない
    			else {
    								// 準備
    				double s  = -A1[p][q];
    				double t  = 0.5 * (A1[p][p] - A1[q][q]);
    				double v  = Math.Abs(t) / Math.Sqrt(s * s + t * t);
    				double sn = Math.Sqrt(0.5 * (1.0 - v));
    				if (s*t < 0.0)
    					sn = -sn;
    				double cs = Math.Sqrt(1.0 - sn * sn);
    								// Akの計算
    				for (int i1 = 0; i1 < n; i1++) {
    					if (i1 == p) {
    						for (int i2 = 0; i2 < n; i2++) {
    							if (i2 == p)
    								A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                               2.0 * A1[p][q] * sn * cs;
    							else if (i2 == q)
    								A2[p][q] = 0.0;
    							else
    								A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
    						}
    					}
    					else if (i1 == q) {
    						for (int i2 = 0; i2 < n; i2++) {
    							if (i2 == q)
    								A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                               2.0 * A1[p][q] * sn * cs;
    							else if (i2 == p)
    								A2[q][p] = 0.0;
    							else
    								A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
    						}
    					}
    					else {
    						for (int i2 = 0; i2 < n; i2++) {
    							if (i2 == p)
    								A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
    							else if (i2 == q)
    								A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
    							else
    								A2[i1][i2] = A1[i1][i2];
    						}
    					}
    				}
    								// Xkの計算
    				for (int i1 = 0; i1 < n; i1++) {
    					for (int i2 = 0; i2 < n; i2++) {
    						if (i2 == p)
    							X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
    						else if (i2 == q)
    							X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
    						else
    							X2[i1][i2] = X1[i1][i2];
    					}
    				}
    								// 次のステップへ
    				k++;
    				for (int i1 = 0; i1 < n; i1++) {
    					for (int i2 = 0; i2 < n; i2++) {
    						A1[i1][i2] = A2[i1][i2];
    						X1[i1][i2] = X2[i1][i2];
    					}
    				}
    			}
    		}
    
    		return ind;
    	}
    }
    
    /*
    ---------データ例(コメント部分を除いて下さい)---------
    4 100   // 変数の数(p)とデータの数(n)
    66 22 44 31   // x1, x2, x3, x4
    25 74 17 81
    50 23 53 71
    25 57 19 81
    74 47 64 47
    39 33 48 46
    14 22 9 69
    67 60 49 26
    42 40 77 65
    11 80 0 86
    32 0 43 74
    68 69 44 68
    24 49 9 71
    42 74 28 46
    60 58 73 28
    36 37 33 68
    24 44 19 83
    30 40 31 50
    55 40 60 49
    63 47 94 41
    72 30 100 45
    19 22 13 75
    43 39 43 34
    90 83 92 31
    51 77 52 82
    53 70 34 31
    28 51 53 44
    40 62 42 79
    31 48 22 68
    57 29 51 30
    64 89 57 42
    49 82 72 29
    53 31 55 43
    79 52 70 10
    45 19 43 57
    35 34 34 89
    4 69 0 100
    49 49 66 66
    92 82 97 6
    5 89 0 100
    65 26 83 28
    56 36 64 38
    48 50 25 22
    30 30 15 55
    40 65 38 42
    14 67 9 67
    84 96 90 8
    53 64 51 54
    50 89 60 52
    76 41 68 9
    49 40 53 49
    78 66 66 17
    76 58 90 29
    41 15 40 49
    63 60 55 33
    40 36 49 67
    78 54 71 18
    62 72 69 12
    64 47 42 53
    56 64 9 15
    77 35 56 25
    44 12 46 87
    80 9 56 19
    36 21 52 78
    48 63 64 48
    43 61 50 47
    58 23 28 50
    90 12 100 0
    13 33 11 77
    67 44 48 28
    75 45 68 17
    81 22 89 9
    46 45 59 55
    56 49 64 55
    65 62 72 27
    34 49 29 77
    45 33 60 63
    20 45 14 99
    33 38 26 87
    44 51 69 52
    64 57 64 48
    44 64 51 28
    63 48 56 11
    29 39 33 84
    40 48 51 54
    40 38 26 62
    68 46 61 26
    58 45 68 48
    64 44 77 63
    59 62 44 66
    81 53 93 19
    23 34 12 68
    51 35 55 46
    74 70 84 17
    42 33 56 44
    46 31 46 53
    33 57 38 63
    40 24 20 42
    53 36 60 31
    0 34 0 100
    */
    			

  8. VB

    '**************************'
    ' 主成分分析               '
    '      coded by Y.Suganuma '
    '**************************'
    Imports System.Text.RegularExpressions
    
    Module Test
    
    	Sub Main()
    
    		Dim MS As Regex = New Regex("\s+") 
    					' 変数の数とデータの数
    		Dim str() As String = MS.Split(Console.ReadLine().Trim())
    		Dim p As Integer    = Integer.Parse(str(0))
    		Dim n As Integer    = Integer.Parse(str(1))
    
    		Dim r(p) As Double
    		Dim x(p,n) As Double
    		Dim a(p,p) As Double
    					' データ
    		For i1 As Integer = 0 To n-1
    			str = MS.Split(Console.ReadLine().Trim())
    			For i2 As Integer = 0 To p-1
    				x(i2,i1) = Double.Parse(str(i2))
    			Next
    		Next
    
    		Dim sw As Integer = principal(p, n, x, r, a, 1.0e-10, 200)
    
    		If sw = 0
    			For i1 As Integer = 0 To p-1
    				Console.Write("主成分 " & r(i1))
    				Console.Write(" 係数")
    				For i2 As Integer = 0 To p-1
    					Console.Write(" " & a(i1,i2))
    				Next
    				Console.WriteLine()
    			Next
    		Else
    			Console.WriteLine("***error  解を求めることができませんでした")
    		End If
    
    	End Sub
    
    	'*********************************'
    	' 主成分分析                      '
    	'      p : 変数の数               '
    	'      n : データの数             '
    	'      x : データ                 '
    	'      r : 分散(主成分)         '
    	'      a : 係数                   '
    	'      eps : 収束性を判定する規準 '
    	'      ct : 最大繰り返し回数      '
    	'      return : =0 : 正常         '
    	'               =1 : エラー       '
    	'*********************************'
    	Function principal(p As Integer, n As Integer, x(,) As Double, r() As Double,
    	                   a(,) As Double, eps As Double, ct As Integer)
    					' 領域の確保
    		Dim C(p,p) As Double
    		Dim A1(p,p) As Double
    		Dim A2(p,p) As Double
    		Dim X1(p,p) As Double
    		Dim X2(p,p) As Double
    					' データの基準化
    		For i1 As Integer = 0 To p-1
    			Dim mean As Double = 0.0
    			Dim s2 As Double   = 0.0
    			For i2 As Integer = 0 To n-1
    				mean += x(i1,i2)
    				s2   += x(i1,i2) * x(i1,i2)
    			Next
    			mean /= n
    			s2   /= n
    			s2    = n * (s2 - mean * mean) / (n - 1)
    			s2    = Math.Sqrt(s2)
    			For i2 As Integer = 0 To n-1
    				x(i1,i2) = (x(i1,i2) - mean) / s2
    			Next
    		Next
    					' 分散強分散行列の計算
    		For i1 As Integer = 0 To p-1
    			For i2 As Integer = i1 To p-1
    				Dim s2 As Double = 0.0
    				For i3 As Integer = 0 To n-1
    					s2 += x(i1,i3) * x(i2,i3)
    				Next
    				s2        /= (n - 1)
    				C(i1,i2)  = s2
    				If i1 <> i2
    					C(i2,i1) = s2
    				End If
    			Next
    		Next
    					' 固有値と固有ベクトルの計算(ヤコビ法)
    		Dim sw As Integer = Jacobi(p, ct, eps, C, A1, A2, X1, X2)
    
    		If sw = 0
    			For i1 As Integer = 0 To p-1
    				r(i1) = A1(i1,i1)
    				For i2 As Integer = 0 To p-1
    					a(i1,i2) = X1(i2,i1)
    				Next
    			Next
    		End If
    
    		Return sw
    
    	End Function
    
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    	' 実対称行列の固有値・固有ベクトル(ヤコビ法)              '
    	'      n : 次数                                             '
    	'      ct : 最大繰り返し回数                                '
    	'      eps : 収束判定条件                                   '
    	'      A : 対象とする行列                                   '
    	'      A1, A2 : 作業域(nxnの行列),A1の対角要素が固有値   '
    	'      X1, X2 : 作業域(nxnの行列),X1の各列が固有ベクトル '
    	'      return : =0 : 正常                                   '
    	'               =1 : 収束せず                               '
    	'      coded by Y.Suganuma                                  '
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    	Function Jacobi(n As Integer, ct As Integer, eps As Double, A(,) As Double,
    	                A1(,) As Double, A2(,) As Double, X1(,) As Double, X2(,) As Double)
    					' 初期設定
    		For i1 As Integer = 0 To n-1
    			For i2 As Integer = 0 To n-1
    				A1(i1,i2) = A(i1,i2)
    				X1(i1,i2) = 0.0
    			Next
    			X1(i1,i1) = 1.0
    		Next
    					' 計算
    		Dim k As Integer   = 0
    		Dim ind As Integer = 1
    		Do While ind > 0 and k < ct
    						' 最大要素の探索
    			Dim max As Double = 0.0
    			Dim p As Integer  = 0
    			Dim q As Integer  = 0
    			For i1 As Integer = 0 To n-1
    				For i2 As Integer = 0 To n-1
    					If i2 <> i1
    						If Math.Abs(A1(i1,i2)) > max
    							max = Math.Abs(A1(i1,i2))
    							p   = i1
    							q   = i2
    						End If
    					End If
    				Next
    			Next
    						' 収束判定
    							' 収束した
    			If max < eps
    				ind = 0
    							' 収束しない
    			Else
    								' 準備
    				Dim s As Double  = -A1(p,q)
    				Dim t As Double  = 0.5 * (A1(p,p) - A1(q,q))
    				Dim v As Double  = Math.Abs(t) / Math.Sqrt(s * s + t * t)
    				Dim sn As Double = Math.Sqrt(0.5 * (1.0 - v))
    				If s*t < 0.0
    					sn = -sn
    				End If
    				Dim cs As Double = Math.Sqrt(1.0 - sn * sn)
    								' Akの計算
    				For i1 As Integer = 0 To n-1
    					If i1 = p
    						For i2 As Integer = 0 To n-1
    							If i2 = p
    								A2(p,p) = A1(p,p) * cs * cs + A1(q,q) * sn * sn -
                                               2.0 * A1(p,q) * sn * cs
    							ElseIf i2 = q
    								A2(p,q) = 0.0
    							Else
    								A2(p,i2) = A1(p,i2) * cs - A1(q,i2) * sn
    							End If
    						Next
    					ElseIf i1 = q
    						For i2 As Integer = 0 To n-1
    							If i2 = q
    								A2(q,q) = A1(p,p) * sn * sn + A1(q,q) * cs * cs +
                                               2.0 * A1(p,q) * sn * cs
    							ElseIf i2 = p
    								A2(q,p) = 0.0
    							Else
    								A2(q,i2) = A1(q,i2) * cs + A1(p,i2) * sn
    							End If
    						Next
    					Else
    						For i2 As Integer = 0 To n-1
    							If i2 = p
    								A2(i1,p) = A1(i1,p) * cs - A1(i1,q) * sn
    							ElseIf i2 = q
    								A2(i1,q) = A1(i1,q) * cs + A1(i1,p) * sn
    							Else
    								A2(i1,i2) = A1(i1,i2)
    							End If
    						Next
    					End If
    				Next
    								' Xkの計算
    				For i1 As Integer = 0 To n-1
    					For i2 As Integer = 0 To n-1
    						If i2 = p
    							X2(i1,p) = X1(i1,p) * cs - X1(i1,q) * sn
    						ElseIf i2 = q
    							X2(i1,q) = X1(i1,q) * cs + X1(i1,p) * sn
    						Else
    							X2(i1,i2) = X1(i1,i2)
    						End If
    					Next
    				Next
    								' 次のステップへ
    				k += 1
    				For i1 As Integer = 0 To n-1
    					For i2 As Integer = 0 To n-1
    						A1(i1,i2) = A2(i1,i2)
    						X1(i1,i2) = X2(i1,i2)
    					Next
    				Next
    			End If
    		Loop
    
    		Return ind
    
    	End Function
    
    End Module
    
    /*
    ---------データ例(コメント部分を除いて下さい)---------
    4 100   // 変数の数(p)とデータの数(n)
    66 22 44 31   // x1, x2, x3, x4
    25 74 17 81
    50 23 53 71
    25 57 19 81
    74 47 64 47
    39 33 48 46
    14 22 9 69
    67 60 49 26
    42 40 77 65
    11 80 0 86
    32 0 43 74
    68 69 44 68
    24 49 9 71
    42 74 28 46
    60 58 73 28
    36 37 33 68
    24 44 19 83
    30 40 31 50
    55 40 60 49
    63 47 94 41
    72 30 100 45
    19 22 13 75
    43 39 43 34
    90 83 92 31
    51 77 52 82
    53 70 34 31
    28 51 53 44
    40 62 42 79
    31 48 22 68
    57 29 51 30
    64 89 57 42
    49 82 72 29
    53 31 55 43
    79 52 70 10
    45 19 43 57
    35 34 34 89
    4 69 0 100
    49 49 66 66
    92 82 97 6
    5 89 0 100
    65 26 83 28
    56 36 64 38
    48 50 25 22
    30 30 15 55
    40 65 38 42
    14 67 9 67
    84 96 90 8
    53 64 51 54
    50 89 60 52
    76 41 68 9
    49 40 53 49
    78 66 66 17
    76 58 90 29
    41 15 40 49
    63 60 55 33
    40 36 49 67
    78 54 71 18
    62 72 69 12
    64 47 42 53
    56 64 9 15
    77 35 56 25
    44 12 46 87
    80 9 56 19
    36 21 52 78
    48 63 64 48
    43 61 50 47
    58 23 28 50
    90 12 100 0
    13 33 11 77
    67 44 48 28
    75 45 68 17
    81 22 89 9
    46 45 59 55
    56 49 64 55
    65 62 72 27
    34 49 29 77
    45 33 60 63
    20 45 14 99
    33 38 26 87
    44 51 69 52
    64 57 64 48
    44 64 51 28
    63 48 56 11
    29 39 33 84
    40 48 51 54
    40 38 26 62
    68 46 61 26
    58 45 68 48
    64 44 77 63
    59 62 44 66
    81 53 93 19
    23 34 12 68
    51 35 55 46
    74 70 84 17
    42 33 56 44
    46 31 46 53
    33 57 38 63
    40 24 20 42
    53 36 60 31
    0 34 0 100
    */
    			

静岡理工科大学 菅沼ホーム 目次 索引