/********************************/ /* 重回帰分析 */ /* coded by Y.Suganuma */ /********************************/ import java.io.*; public class Test { public static void main(String args[]) throws IOException { double b[], y[], X[][]; int i1, i2, n, N, sw; String str; Data dt; BufferedReader in = new BufferedReader(new InputStreamReader(System.in)); // 説明変数の数とデータの数 str = in.readLine(); dt = new Data(str); dt.next(); n = Integer.parseInt(dt.data); dt.next(); N = Integer.parseInt(dt.data); y = new double [N]; X = new double [N][n+1]; b = new double [n+1]; for (i1 = 0; i1 < N; i1++) { // データ X[i1][0] = 1.0; str = in.readLine(); dt = new Data(str); dt.next(); y[i1] = Double.parseDouble(dt.data); for (i2 = 0; i2 < n; i2++) { dt.next(); X[i1][i2+1] = Double.parseDouble(dt.data); } } sw = App.regression(n, N, X, y, b, 1.0e-10); if (sw == 0) { System.out.println("結果"); for (i1 = 0; i1 < n+1; i1++) System.out.println(" b" + i1 + " " + b[i1]); } else System.out.println("***error 逆行列を求めることができませんでした"); } } /**************************** /* 科学技術系算用の手法 */ /****************************/ class App { /***************************************/ /* 重回帰分析 */ /* n : 説明変数の数 */ /* N : データの数 */ /* X,y : データ */ /* b : 偏回帰係数 */ /* eps : 正則性を判定する規準 */ /* return : =0 : 正常 */ /* =1 : エラー */ /***************************************/ static int regression(int n, int N, double X[][], double y[], double b[], double eps) { double w[][]; int i1, i2, i3, sw; n++; w = new double [n][n+1]; for (i1 = 0; i1 < n; i1++) { for (i2 = 0; i2 < n; i2++) { w[i1][i2] = 0.0; for (i3 = 0; i3 < N; i3++) w[i1][i2] += X[i3][i1] * X[i3][i2]; } } for (i1 = 0; i1 < n; i1++) { w[i1][n] = 0.0; for (i2 = 0; i2 < N; i2++) w[i1][n] += X[i2][i1] * y[i2]; } sw = App.gauss(w, n, 1, eps); if (sw == 0) { for (i1 = 0; i1 < n; i1++) b[i1] = w[i1][n]; } else sw = 1; return sw; } /***********************************************************/ /* 線形連立方程式を解く(逆行列を求める) */ /* 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------------------------- /**********************************************/ /* 1行のデータから個々のデータを取り出す */ /* 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; } } -------------------------データ------------------------- 3 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