/****************************/
/* 重回帰分析 */
/* coded by Y.Suganuma */
/****************************/
#include <stdio.h>
#include <math.h>
int gauss(double **, int, int, double);
double *regression(int, int, double **, double *, double);
int main()
{
double *b, *y, **X;
int i1, i2, n, N;
scanf("%d %d", &n, &N); // 説明変数の数とデータの数
y = new double [N];
X = new double * [N];
for (i1 = 0; i1 < N; i1++)
X[i1] = new double [n+1];
for (i1 = 0; i1 < N; i1++) { // データ
X[i1][0] = 1.0;
scanf("%lf", &y[i1]);
for (i2 = 0; i2 < n; i2++)
scanf("%lf", &X[i1][i2+1]);
}
b = regression(n, N, X, y, 1.0e-10);
if (b != NULL) {
printf("結果\n");
for (i1 = 0; i1 < n+1; i1++)
printf(" b%d %f\n", i1, b[i1]);
delete [] b;
}
else
printf("***error 逆行列を求めることができませんでした\n");
for (i1 = 0; i1 < N; i1++)
delete [] X[i1];
delete [] X;
delete [] y;
return 0;
}
/******************************************/
/* 重回帰分析 */
/* n : 説明変数の数 */
/* N : データの数 */
/* X,y : データ */
/* eps : 正則性を判定する規準 */
/* return : 偏回帰係数 */
/* エラーの場合はNULLを返す */
/******************************************/
double *regression(int n, int N, double **X, double *y, double eps)
{
double **w, *b;
int i1, i2, i3, sw;
n++;
b = new double [n];
w = new double * [n];
for (i1 = 0; i1 < n; i1++)
w[i1] = new double [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 = gauss(w, n, 1, eps);
if (sw == 0) {
for (i1 = 0; i1 < n; i1++)
b[i1] = w[i1][n];
}
else
b = NULL;
for (i1 = 0; i1 < n; i1++)
delete [] w[i1];
delete [] w;
return b;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
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 = fabs(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);
}
---------データ例（コメント部分を除いて下さい）---------
3 100 // 説明変数の数（n）とデータの数（N）
66 22 44 31 // y, x1, x2, x3
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