#概要
####最小二乗法
最小二乗法は、誤差を伴う測定値の処理において、その誤差の二乗の和を最小にすることで、関係式を求める方法です。
####プログラムについて
今回作成したプログラムは、予め用意されたデータx
y
をリストにし、そのデータから行列S
T
を求め、係数a
を導く処理となります。また、求めたい多項式の次数をm
と置いています。
最小二乗法を簡単に求めることができるライブラリは既に存在するかもしれませんが、今回はそのようなライブラリは利用せず、配列操作ライブラリnumpy
のみを利用し、データ量を変更しても融通が利くような処理を作成しました。
#プログラム
# -*- coding: utf-8 -*-
#!/usr/bin/env python3
import numpy as np
X_data = [10, 20, 50, 70, 100]
Y_data = [9.3, 9.8, 10.9, 11.9, 13.1]
m = 1
def ST_list(power):
S, T = 0, 0
for i, j in zip(X_data,Y_data):
S += i**power
T += i**power * j
return S, T
def A_matrix(n_size, len):
A = []
for k in range(n_size):
A.append(S_list[k+len-m*2-1:k+len-m*2-1+n_size])
return A
S_list = []
T_list = []
for n in range(len(X_data)):
S , T = ST_list(n)
S_list.insert(0, S)
T_list.insert(0, T)
A = np.array(A_matrix(m+1, len(S_list)))
b = np.array(T_list[-1*m-1::]).reshape(m+1, 1)
Ainv = np.linalg.inv(A)
x = np.dot(Ainv, b)
print("S ="); print(A); print()
print("T ="); print(b); print()
print("a ="); print(x)
#実行結果
#####一次多項式の係数aを求める
m = 1
S =
[[17900 250]
[ 250 5]]
T =
[[2977.]
[ 55.]]
a =
[[0.04203704]
[8.89814815]]
#####二次多項式の係数aを求める
m = 2
S =
[[130430000 1477000 17900]
[ 1477000 17900 250]
[ 17900 250 5]]
T =
[[2.2141e+05]
[2.9770e+03]
[5.5000e+01]]
a =
[[1.22729504e-05]
[4.07142857e-02]
[8.92034855e+00]]
#まとめ
今回のデータでは、一次多項式と二次多項式のときの係数a
を求めました。データ量を増やせば三次、四次も行うことができますが、二次多項式のように浮動小数の影響を受けるため、見やすい値にするにはround
やformat
を利用する必要がありそうです。
#終わりに
今回の記事はここまでになります。
見にくいプログラムで申し訳ないです。特に配列操作に関してはもっと良い方法があるように感じます。アドバイス等あれば、コメント欄にてお願いいたします。
最後までお付き合い頂きありがとうございました。