More than 5 years have passed since last update.

Numpy.polyfit を使ったカーブフィッティング

Last updated at 2019-05-10Posted at 2019-05-10

様々な補間法と最小2乗法をPythonで理解するのうち、「Numpy.polyfit を使ったカーブフィッティング」を、実データっぽい模擬データを解析するように書き直したサンプルプログラムです。

実データっぽい模擬データ

こういうデータを多項式近似したいとしましょう。

x_observed = [9, 28, 38, 58, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198, 208, 218, 228, 238, 278, 288, 298]
y_observed = [51, 80, 112, 294, 286, 110, 59, 70, 56, 70, 104, 59, 59, 72, 87, 99, 64, 60, 74, 151, 157, 57, 83]

まずはデータの図示

こういうデータですね。

%matplotlib inline
import matplotlib.pyplot as plt

plt.scatter(x_observed, y_observed)
plt.grid()

フィッティンブカーブを描く x 座標を作る

実データの最小値から最大値までの間に等間隔に１００ポイント作りました。

import numpy as np
x_latent = np.linspace(min(x_observed), max(x_observed), 100)

Numpy.polyfit を用いた最小二乗法

たとえば１次式から９次式までフィッティングしてみましょう。

cf1 = ["最小2乗法（1次式）", lambda x, y: np.polyfit(x, y, 1)]
cf2 = ["最小2乗法（2次式）", lambda x, y: np.polyfit(x, y, 2)]
cf3 = ["最小2乗法（3次式）", lambda x, y: np.polyfit(x, y, 3)]
cf4 = ["最小2乗法（4次式）", lambda x, y: np.polyfit(x, y, 4)]
cf5 = ["最小2乗法（5次式）", lambda x, y: np.polyfit(x, y, 5)]
cf6 = ["最小2乗法（6次式）", lambda x, y: np.polyfit(x, y, 6)]
cf7 = ["最小2乗法（7次式）", lambda x, y: np.polyfit(x, y, 7)]
cf8 = ["最小2乗法（8次式）", lambda x, y: np.polyfit(x, y, 8)]
cf9 = ["最小2乗法（9次式）", lambda x, y: np.polyfit(x, y, 9)]

Sympy を用いた数式の表示の準備

Sympy の使い方は微分や微分方程式をPythonで理解するなどを参考にしていただければと。


import sympy as sym
from sympy.plotting import plot
sym.init_printing(use_unicode=True)

x, y = sym.symbols("x y")

最小二乗法によるカーブフィッティング（１次式から９次式まで）

数式の表示とグラフの表示を一気に行います。

for method_name, method in [cf1, cf2, cf3, cf4, cf5, cf6, cf7, cf8, cf9]:
    print(method_name)
    # 係数の計算
    coefficients = method(x_observed, y_observed)
    
    # Sympy を用いた数式の表示
    expr = 0
    for index, coefficient in enumerate(coefficients):
        expr += coefficient * x ** (len(coefficients) - index - 1)
    display(sym.Eq(y, expr))
    
    # プロットと曲線の表示
    fitted_curve = np.poly1d(method(x_observed, y_observed))(x_latent)
    plt.scatter(x_observed, y_observed, label="observed")
    plt.plot(x_latent, fitted_curve, c="red", label="fitted")
    plt.grid()
    plt.legend()
    plt.show()