Help us understand the problem. What is going on with this article?

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

概要

最小二乗法

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

プログラムについて

今回作成したプログラムは、予め用意されたデータ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を求めました。データ量を増やせば三次、四次も行うことができますが、二次多項式のように浮動小数の影響を受けるため、見やすい値にするにはroundformatを利用する必要がありそうです。

終わりに

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

AlexG00t
大学生です。備忘録としてpythonを中心に自分で作成したプログラムを載せていきます。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした