最小二乗法のアルゴリズムについて

前提・実現したいこと

座標データをテキストファイルで読み込み
最小二乗法で近似曲線を計算プログラムを書いています。

発生している問題・エラーメッセージ

アルゴリズムが間違っているみたいで挙動がおかしいです。

該当のソースコード

下記に示したのが最小二乗法の計算アルゴリズムです。

C
1//ガウス消去法（連立方程式を解くため）
2int gauss(double **w, int n, int m, double eps){
3
4	double y1;
5  double y2;
6	int ind = 0;
7  int nm;
8  int m1;
9  int m2;
10  int i;
11  int j;
12  int k;
13
14	nm = n + m;
15
16	for (i=0;i<n&&ind==0;i++) {
17
18		y1 = 0;
19		m1 = i+1;
20		m2 = 0;
21
22		for (j=i;j<n;j++) {
23			y2 = fabs(w[j][i]);
24			if (y1<y2) {
25				y1 = y2;
26				m2 = j;
27			}
28		}
29
30		if (y1<eps){
31			ind = 1;
32		}else{
33			for (j=i;j<nm;j++) {
34				y1 = w[i][j];
35				w[i][j] = w[m2][j];
36				w[m2][j] = y1;
37			}
38
39			y1 = 1.0 / w[i][i];
40
41			for (j=m1;j<nm;j++){
42				w[i][j] *= y1;
43      }
44
45			for (j=0;j<n;j++) {
46				if (j!=i) {
47					for (k=m1;k<nm;k++){
48						w[j][k] -= w[j][i] * w[i][k];
49          }
50				}
51			}
52		}
53	}
54
55	return ind;
56}
57
58//最小二乗法の関数
59double *least_squares_method(int m,int n,double *x,double *y){
60  double **a;//可変長の配列
61  double **w;//可変長の配列
62  double *z;//可変長の配列
63  double x1;
64  double x2;
65
66  int sw;
67  int i;
68  int j;
69  int k;
70
71  m++;
72  z = new double[m];
73  w = new double*[m];
74
75  for (i=0;i<m;i++){
76  w[i] = new double [m+1];
77  }
78  a = new double * [n];
79
80  for (i=0;i<n;i++) {
81		a[i] = new double [m];
82		a[i][m-2] = x[i];
83		a[i][m-1] = 1.0;
84		x1 = a[i][m-2];
85		x2 = x1;
86		for (j=m-3;j>=0;j--){
87			x2 *= x1;
88			a[i][j] = x2;
89		}
90	}
91
92  for (i=0;i<m;i++){
93		for (j=0;j<m;j++) {
94			w[i][j] = 0.0;
95			for (k=0;k<n;k++){
96				w[i][j] += a[k][i]*a[k][j];
97      }
98		}
99	}
100
101  for (i=0;i<m;i++) {
102		w[i][m] = 0.0;
103		for(int j=0;j<n;j++){
104			w[i][m] += a[j][i] * y[j];
105	  }
106  }
107
108  sw=gauss(w,m,1,1.0e-10);
109
110	if (sw==0) {
111		for (i=0;i<m;i++){
112			z[i] = w[i][m];
113    }
114	}else{
115		z = NULL;
116  }
117
118  return z;
119}
120