/********************************/ /* 正準相関分析 */ /* coded by Y.Suganuma */ /********************************/ import java.io.*; class Test { static void main(String args[]) throws IOException { double x[][], ryz[], ab[][], v0[]; int i1, i2, q, r, s, n, sw; String str; Data dt; BufferedReader in = new BufferedReader(new InputStreamReader(System.in)); // 各組の変数の数とデータの数 str = in.readLine(); dt = new Data(str); dt.next(); r = Integer.parseInt(dt.data); dt.next(); s = Integer.parseInt(dt.data); dt.next(); n = Integer.parseInt(dt.data); q = r + s; ryz = new double [r]; v0 = new double [r]; x = new double [q][n]; ab = new double [r][q]; for (i1 = 0; i1 < r; i1++) v0[i1] = 1.0; // データ for (i1 = 0; i1 < n; i1++) { str = in.readLine(); dt = new Data(str); for (i2 = 0; i2 < q; i2++) { dt.next(); x[i2][i1] = Double.parseDouble(dt.data); } } sw = App.canonical(r, s, n, x, ryz, ab, 1.0e-10, v0, 15, 200); if (sw > 0) { for (i1 = 0; i1 < sw; i1++) { System.out.println("相関係数 " + ryz[i1]); System.out.print(" a"); for (i2 = 0; i2 < r; i2++) System.out.print(" " + ab[i1][i2]); System.out.println(); System.out.print(" b"); for (i2 = 0; i2 < s; i2++) System.out.print(" " + ab[i1][r+i2]); System.out.println(); } } else System.out.println("***error 解を求めることができませんでした"); } } /**************************** /* 科学技術系算用の手法 */ /****************************/ class App { /***************************************/ /* 正準相関分析 */ /* r,s : 各組における変数の数 */ /* n : データの数 */ /* x : データ */ /* ryz : 相関係数 */ /* ab : 各組の係数（a,bの順） */ /* eps : 正則性を判定する規準 */ /* v0 : 固有ベクトルの初期値 */ /* m : 丸め誤差の修正回数 */ /* ct : 最大繰り返し回数 */ /* return : >0 : 相関係数の数 */ /* =0 : エラー */ /***************************************/ static int canonical(int r, int s, int n, double x[][], double ryz[], double ab[][], double eps, double v0[], int m, int ct) { double A[][], B[][], C[][], C11[][], C11i[][], C12[][], C21[][], C22[][], C22i[][], mean[], vv, w[][], w1[], len, v1[], v2[]; int i1, i2, i3, q, sw = 0, sw1; // 領域の確保 q = r + s; mean = new double [q]; w1 = new double [s]; v1 = new double [r]; v2 = new double [r]; A = new double [s][s]; B = new double [r][r]; C = new double [q][q]; C11 = new double [r][r]; C11i = new double [r][r]; C12 = new double [r][s]; C21 = new double [s][r]; C22 = new double [s][s]; C22i = new double [s][s]; w = new double [s][2*s]; // 平均値の計算 for (i1 = 0; i1 < q; i1++) { mean[i1] = 0.0; for (i2 = 0; i2 < n; i2++) mean[i1] += x[i1][i2]; mean[i1] /= n; } // 分散強分散行列の計算 for (i1 = 0; i1 < q; i1++) { for (i2 = i1; i2 < q; i2++) { vv = 0.0; for (i3 = 0; i3 < n; i3++) vv += (x[i1][i3] - mean[i1]) * (x[i2][i3] - mean[i2]); vv /= (n - 1); C[i1][i2] = vv; if (i1 != i2) C[i2][i1] = vv; } } // C11, C12, C21, C22 の設定 // C12 for (i1 = 0; i1 < r; i1++) { for (i2 = 0; i2 < s; i2++) C12[i1][i2] = C[i1][i2+r]; } // C21 for (i1 = 0; i1 < s; i1++) { for (i2 = 0; i2 < r; i2++) C21[i1][i2] = C[i1+r][i2]; } // C11とその逆行列 for (i1 = 0; i1 < r; i1++) { for (i2 = 0; i2 < r; i2++) { w[i1][i2] = C[i1][i2]; C11[i1][i2] = C[i1][i2]; } for (i2 = r; i2 < 2*r; i2++) w[i1][i2] = 0.0; w[i1][i1+r] = 1.0; } sw1 = App.gauss(w, r, r, 1.0e-10); if (sw1 == 0) { for (i1 = 0; i1 < r; i1++) { for (i2 = 0; i2 < r; i2++) C11i[i1][i2] = w[i1][i2+r]; } } else sw = 1; // C22とその逆行列 for (i1 = 0; i1 < s; i1++) { for (i2 = 0; i2 < s; i2++) { w[i1][i2] = C[i1+r][i2+r]; C22[i1][i2] = C[i1+r][i2+r]; } for (i2 = s; i2 < 2*s; i2++) w[i1][i2] = 0.0; w[i1][i1+s] = 1.0; } sw1 = App.gauss(w, s, s, eps); if (sw1 == 0) { for (i1 = 0; i1 < s; i1++) { for (i2 = 0; i2 < s; i2++) C22i[i1][i2] = w[i1][i2+s]; } } else sw = 1; // 固有値λ及び固有ベクトルaの計算 if (sw == 0) { // 行列の計算 for (i1 = 0; i1 < s; i1++) { for (i2 = 0; i2 < r; i2++) { A[i1][i2] = 0.0; for (i3 = 0; i3 < s; i3++) A[i1][i2] += C22i[i1][i3] * C21[i3][i2]; } } for (i1 = 0; i1 < r; i1++) { for (i2 = 0; i2 < s; i2++) { w[i1][i2] = 0.0; for (i3 = 0; i3 < s; i3++) w[i1][i2] += C12[i1][i3] * A[i3][i2]; } } for (i1 = 0; i1 < r; i1++) { for (i2 = 0; i2 < r; i2++) { A[i1][i2] = 0.0; for (i3 = 0; i3 < r; i3++) A[i1][i2] += C11i[i1][i3] * w[i3][i2]; } } // 固有値と固有ベクトル（べき乗法） for (i1 = 0; i1 < r; i1++) { for (i2 = 0; i2 < r; i2++) B[i1][i2] = 0.0; B[i1][i1] = 1.0; } sw = App.power(0, r, m, ct, eps, A, B, C, ryz, ab, v0, v1, v2); if (sw > 0) { for (i1 = 0; i1 < r; i1++) { // 相関係数 ryz[i1] = Math.sqrt(ryz[i1]); // 大きさの調整(a) for (i2 = 0; i2 < r; i2++) { w1[i2] = 0.0; for (i3 = 0; i3 < r; i3++) w1[i2] += C11[i2][i3] * ab[i1][i3]; } len = 0.0; for (i2 = 0; i2 < r; i2++) len += ab[i1][i2] * w1[i2]; len = Math.sqrt(len); for (i2 = 0; i2 < r; i2++) ab[i1][i2] /= len; // bの計算 for (i2 = 0; i2 < s; i2++) { w1[i2] = 0.0; for (i3 = 0; i3 < r; i3++) w1[i2] += C21[i2][i3] * ab[i1][i3]; } for (i2 = 0; i2 < s; i2++) { ab[i1][i2+r] = 0.0; for (i3 = 0; i3 < s; i3++) ab[i1][i2+r] += C22i[i2][i3] * w1[i3]; } for (i2 = 0; i2 < s; i2++) ab[i1][i2+r] /= ryz[i1]; // 大きさの調整(b) for (i2 = 0; i2 < s; i2++) { w1[i2] = 0.0; for (i3 = 0; i3 < s; i3++) w1[i2] += C22[i2][i3] * ab[i1][i3+r]; } len = 0.0; for (i2 = 0; i2 < s; i2++) len += ab[i1][i2+r] * w1[i2]; len = Math.sqrt(len); for (i2 = 0; i2 < s; i2++) ab[i1][i2+r] /= len; } } } else sw = 0; return sw; } /********************************************/ /* 最大固有値と固有ベクトル（べき乗法） */ /* i : 何番目の固有値かを示す */ /* n : 次数 */ /* m : 丸め誤差の修正回数 */ /* ct : 最大繰り返し回数 */ /* eps : 収束判定条件 */ /* A : 対象とする行列 */ /* B : nxnの行列，最初は，単位行列 */ /* C : 作業域，nxnの行列 */ /* r : 固有値 */ /* v : 各行が固有ベクトル（nxn） */ /* v0 : 固有ベクトルの初期値 */ /* v1,v2 : 作業域（n次元ベクトル） */ /* return : 求まった固有値の数 */ /* coded by Y.Suganuma */ /********************************************/ static int power(int i, int n, int m, int ct, double eps, double A[][], double B[][], double C[][], double r[], double v[][], double v0[], double v1[], double v2[]) { double l1, l2, x; int i1, i2, i3, ind, k, k1; // 初期設定 ind = i; k = 0; l1 = 0.0; for (i1 = 0; i1 < n; i1++) { l1 += v0[i1] * v0[i1]; v1[i1] = v0[i1]; } l1 = Math.sqrt(l1); // 繰り返し計算 while (k < ct) { // 丸め誤差の修正 if (k%m == 0) { l2 = 0.0; for (i1 = 0; i1 < n; i1++) { v2[i1] = 0.0; for (i2 = 0; i2 < n; i2++) v2[i1] += B[i1][i2] * v1[i2]; l2 += v2[i1] * v2[i1]; } l2 = Math.sqrt(l2); for (i1 = 0; i1 < n; i1++) v1[i1] = v2[i1] / l2; } // 次の近似 l2 = 0.0; for (i1 = 0; i1 < n; i1++) { v2[i1] = 0.0; for (i2 = 0; i2 < n; i2++) v2[i1] += A[i1][i2] * v1[i2]; l2 += v2[i1] * v2[i1]; } l2 = Math.sqrt(l2); for (i1 = 0; i1 < n; i1++) v2[i1] /= l2; // 収束判定 // 収束した場合 if (Math.abs((l2-l1)/l1) < eps) { k1 = -1; for (i1 = 0; i1 < n && k1 < 0; i1++) { if (Math.abs(v2[i1]) > 0.001) { k1 = i1; if (v2[k1]*v1[k1] < 0.0) l2 = -l2; } } k = ct; r[i] = l2; for (i1 = 0; i1 < n; i1++) v[i][i1] = v2[i1]; if (i == n-1) ind = i + 1; else { for (i1 = 0; i1 < n; i1++) { for (i2 = 0; i2 < n; i2++) { C[i1][i2] = 0.0; for (i3 = 0; i3 < n; i3++) { x = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3]; C[i1][i2] += x * B[i3][i2]; } } } for (i1 = 0; i1 < n; i1++) { for (i2 = 0; i2 < n; i2++) B[i1][i2] = C[i1][i2]; } for (i1 = 0; i1 < n; i1++) { v1[i1] = 0.0; for (i2 = 0; i2 < n; i2++) v1[i1] += B[i1][i2] * v0[i2]; } for (i1 = 0; i1 < n; i1++) v0[i1] = v1[i1]; ind = App.power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2); } } // 収束しない場合 else { for (i1 = 0; i1 < n; i1++) v1[i1] = v2[i1]; l1 = l2; k++; } } return ind; } /***********************************************************/ /* 線形連立方程式を解く（逆行列を求める） */ /* w : 方程式の左辺及び右辺 */ /* n : 方程式の数 */ /* m : 方程式の右辺の列の数 */ /* eps : 逆行列の存在を判定する規準 */ /* return : =0 : 正常 */ /* =1 : 逆行列が存在しない */ /***********************************************************/ static int gauss(double w[][], int n, int m, double eps) { double y1, y2; int ind = 0, nm, m1, m2, i1, i2, i3; nm = n + m; for (i1 = 0; i1 < n && ind == 0; i1++) { y1 = .0; m1 = i1 + 1; m2 = 0; // ピボット要素の選択 for (i2 = i1; i2 < n; i2++) { y2 = Math.abs(w[i2][i1]); if (y1 < y2) { y1 = y2; m2 = i2; } } // 逆行列が存在しない if (y1 < eps) ind = 1; // 逆行列が存在する else { // 行の入れ替え for (i2 = i1; i2 < nm; i2++) { y1 = w[i1][i2]; w[i1][i2] = w[m2][i2]; w[m2][i2] = y1; } // 掃き出し操作 y1 = 1.0 / w[i1][i1]; for (i2 = m1; i2 < nm; i2++) w[i1][i2] *= y1; for (i2 = 0; i2 < n; i2++) { if (i2 != i1) { for (i3 = m1; i3 < nm; i3++) w[i2][i3] -= w[i2][i1] * w[i1][i3]; } } } } return(ind); } } -------------------------class Data------------------------- /**********************************************/ /* １行のデータから個々のデータを取り出す */ /* coded by Y.Suganuma */ /**********************************************/ public class Data { private int len, start, end; private String line; String data; // コンストラクタ Data (String l) { line = l; len = line.length(); start = 0; end = 0; data = new String (); // 次のデータ } // 次のデータの取り出し // =0 : 成功，=-1 : 失敗 int next() { int k = 0; if (start < len) { int sw = 0; while (sw == 0) { end = line.indexOf(" ", start); if (end >= 0) { if ((end - start) > 0) { data = line.substring(start, end); sw = 1; } } else { end = len; sw = 1; if ((end - start) > 0) data = line.substring(start, end); else k = -1; } start = end + 1; } } else k = -1; return k; } } -------------------------データ------------------------- 2 2 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