今回は、線形一次方程式のカーブフィッティングを行う方法を紹介します。

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]