Python NumPy サンプルコード: 線形一次方程式のカーブフィッティング
今回は、線形一次方程式のカーブフィッティングを行う方法を紹介します。
Python で線形一次方程式のカーブフィッティング
線形一次方程式のカーブフィッティングは最小二乗法による polyfit 関数を使うのが簡単です。 しかし、機械工学に関していうならば、動特性の同定やロボットの制御などで良く用いられるムーア-ペンローズの擬似逆行列による方法も知っておくといいです。
線形一次方程式を
\begin{equation} ax+b = y \nonumber \end{equation}とし、n 個の測定データがこれによく当てはまるとすると、以下の行列形式の方程式が得られます。
\begin{equation} \begin{bmatrix} x_1 & 1\\ \vdots & \vdots\\ x_n & 1 \end{bmatrix}\begin{bmatrix} {a}\\ {b} \end{bmatrix} = \begin{bmatrix} y_1\\ \vdots\\ y_n \end{bmatrix} \nonumber \end{equation}ムーア-ペンローズの擬似逆行列を用いれば、パラメータ \(a, b\) の最適値を推定できます。
\begin{equation} \begin{bmatrix} {a}\\ {b} \end{bmatrix} = \begin{bmatrix} x_1 & 1\\ \vdots & \vdots\\ x_n & 1 \end{bmatrix}^+\begin{bmatrix} y_1\\ \vdots\\ y_n \end{bmatrix} \nonumber \end{equation}Python では以下のように線形一次方程式でのカーブフィッティング問題を解くことができます。
# coding: utf-8 import numpy as np import matplotlib.pyplot as plt # polyfit x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) y = np.array([0.2, 0.3, 0.7, 0.6, 0.8, 1.0]) p1 = np.polyfit(x, y, 1) print("p1=", p1) # pinv x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) y = np.array([0.2, 0.3, 0.7, 0.6, 0.8, 1.0]) A = np.c_[x, np.ones_like(x)] p2 = np.dot(np.linalg.pinv(A), y) print("p2=", p2) # lstsq x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) y = np.array([0.2, 0.3, 0.7, 0.6, 0.8, 1.0]) A = np.c_[x, np.ones_like(x)] p3 = np.linalg.lstsq(A, y)[0] print("p3=", p3) xp = np.linspace(0, 5, 100) plt.figure() plt.plot(x, y, 'r.') plt.plot(xp, p2[0]*xp+p2[1], 'b') # plt.show()
p1= [ 0.15428571 0.21428571] p2= [ 0.15428571 0.21428571] p3= [ 0.15428571 0.21428571]