概要

最小二乗法

最小二乗法は、誤差を伴う測定値の処理において、その誤差の二乗の和を最小にすることで、関係式を求める方法です。

プログラムについて

今回作成したプログラムは、予め用意されたデータx yをリストにし、そのデータから行列S Tを求め、係数aを導く処理となります。また、求めたい多項式の次数をmと置いています。
最小二乗法を簡単に求めることができるライブラリは既に存在するかもしれませんが、今回はそのようなライブラリは利用せず、配列操作ライブラリnumpyのみを利用し、データ量を変更しても融通が利くような処理を作成しました。

プログラム

LeastSquares.py

# -*- 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を利用する必要がありそうです。

終わりに

今回の記事はここまでになります。
見にくいプログラムで申し訳ないです。特に配列操作に関してはもっと良い方法があるように感じます。アドバイス等あれば、コメント欄にてお願いいたします。
最後までお付き合い頂きありがとうございました。

最小二乗法で多項式の係数を求める

概要