1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
| /*********************************************
* 最小二乗法 *
*********************************************/
#include <iostream> // for cout
#include <stdio.h> // for printf()
#define N 7 // データ数
#define M 5 // 予測曲線の次数
using namespace std;
/*
* 計算クラス
*/
class Calc
{
// 各種変数
double a[M + 1][M + 2], s[2 * M + 1], t[M + 1];
int i, j, k;
double p, d, px;
public:
Calc(); // コンストラクタ
void calcLeastSquaresMethod(); // 最小二乗法
double ipow(double, int); // べき乗計算
};
/*
* コンストラクタ
*/
Calc::Calc()
{
// s, t 初期化
for (i = 0; i <= 2 * M; i++)
s[i] = 0;
for (i = 0; i <= M; i++)
t[i] = 0;
}
/*
* 最小二乗法
*/
void Calc::calcLeastSquaresMethod()
{
// 測定データ
static double x[] = {-3, -2, -1, 0, 1, 2, 3};
static double y[] = { 5, -2, -3, -1, 1, 4, 5};
// s[], t[] 計算
for (i = 0; i < N; i++) {
for (j = 0; j <= 2 * M; j++)
s[j] += ipow(x[i], j);
for (j = 0; j <= M; j++)
t[j] += ipow(x[i], j) * y[i];
}
// a[][] に s[], t[] 代入
for (i = 0; i <= M; i++) {
for (j = 0; j <= M; j++)
a[i][j] = s[i + j];
a[i][M + 1] = t[i];
}
// 掃き出し
for (k = 0; k <= M; k++) {
p = a[k][k];
for (j = k; j <= M + 1; j++)
a[k][j] /= p;
for (i = 0; i <= M; i++) {
if (i != k) {
d = a[i][k];
for (j = k; j <= M + 1; j++)
a[i][j] -= d * a[k][j];
}
}
}
// y 値計算＆結果出力
for (k = 0; k <= M; k++)
printf("a%d = %10.6f\n", k, a[k][M + 1]);
printf(" x y\n");
for (px = -3; px <= 3; px += .5) {
p = 0;
for (k = 0; k <= M; k++)
p += a[k][M + 1] * ipow(px, k);
printf("%5.1f%5.1f\n", px, p);
}
}
/*
* べき乗計算
*/
double Calc::ipow(double p, int n)
{
int k;
double s = 1;
for (k = 1; k <= n; k++)
s *= p;
return s;
}
/*
* メイン処理
*/
int main()
{
try
{
// 計算クラスインスタンス化
Calc objCalc;
// 最小二乗法
objCalc.calcLeastSquaresMethod();
}
catch (...) {
cout << "例外発生！" << endl;
return -1;
}
// 正常終了
return 0;
}
|
