Python - 最小二乗法！
Updated:
今回は、近似方程式を「最小二乗法」で解くアルゴリズムを Python3 で実装してみました。
0. 前提条件
- LMDE 2 (Linux Mint Debian Edition 2; 64bit) での作業を想定。
- Python 3.6.4 での作業を想定。
- 当方は他のバージョンとの共存環境であり、
python3.6,pip3.6で 3.6 系を使用するようにしている。（適宜、置き換えて考えること）
1. アルゴリズムについて
当ブログ過去記事を参照。
2. Python スクリプトの作成
- 敢えてオブジェクト指向で作成している。
- Shebang ストリング（1行目）では、フルパスでコマンド指定している。（当方の慣習）
- 数値計算ライブラリ NumPy は使用しない。（この程度の行列計算は List で充分）
- 必要であれば、スクリプト内の定数を変更する。（解きたい近似方程式に合わせて）
File: least_squares_method.py
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
#! /usr/local/bin/python3.6
"""
Solving approximate equation with least squares method
"""
import sys
import traceback
class LeastSquaresMethod:
N, M = 7, 5 # Number of data, Dimension of curve
X = [-3, -2, -1, 0, 1, 2, 3] # Measurement data x
Y = [ 5, -2, -3, -1, 1, 4, 5] # Measurement data y
def __init__(self):
self.a = [[0.0 for _ in range(self.M + 2)] for _ in range(self.M + 1)]
self.s = [0.0 for _ in range(2 * self.M + 1)]
self.t = [0.0 for _ in range(self.M + 1)]
def exec(self):
""" Executeion """
try:
self.__calc_st()
self.__ins_st()
self.__sweep_out()
self.__display()
except Exception as e:
raise
def __calc_st(self):
""" Calculation of s and t """
try:
for i in range(self.N):
for j in range(2 * self.M + 1):
self.s[j] += (self.X[i] ** j)
for j in range(self.M + 1):
self.t[j] += (self.X[i] ** j) * self.Y[i]
except Exception as e:
raise
def __ins_st(self):
""" Insertion of s and t into a """
try:
for i in range(self.M + 1):
for j in range(self.M + 1):
self.a[i][j] = self.s[i + j]
self.a[i][self.M + 1] = self.t[i]
except Exception as e:
raise
def __sweep_out(self):
""" Sweeping out """
try:
for k in range(self.M + 1):
p = self.a[k][k]
for j in range(k, self.M + 2):
self.a[k][j] /= p
for i in range(self.M + 1):
if i == k:
continue
d = self.a[i][k]
for j in range(k, self.M + 2):
self.a[i][j] -= d * self.a[k][j]
except Exception as e:
raise
def __display(self):
""" Computation of y value and display """
try:
for k in range(self.M + 1):
print("a{:d} = {:10.6f}".format(k, self.a[k][self.M + 1]))
print(" x y")
for px in [0.5 * i for i in range(-6, 7)]:
p = 0
for k in range(self.M + 1):
p += self.a[k][self.M + 1] * (px ** k)
print("{:5.1f}{:5.1f}".format(px, p))
except Exception as e:
raise
if __name__ == '__main__':
try:
obj = LeastSquaresMethod()
obj.exec()
except Exception as e:
traceback.print_exc()
sys.exit(1)
3. Python スクリプトの実行
まず、実行権限を付与。
$ chmod +x least_squares_method.py
そして、実行。
$ ./least_squares_method.py
a0 = -1.259740
a1 = 2.100000
a2 = 0.424242
a3 = -0.083333
a4 = 0.030303
a5 = -0.016667
x y
-3.0 5.0
-2.5 0.3
-2.0 -2.1
-1.5 -2.9
-1.0 -2.8
-0.5 -2.2
0.0 -1.3
0.5 -0.1
1.0 1.2
1.5 2.6
2.0 3.9
2.5 4.9
3.0 5.0
C++ 版等と同じ結果になるはず。
以上
Comments