【Python/Pandas】信号データの近似直線・曲線を計算【Scikit-learn】

PythonモジュールPandasで読み込んだ信号データの近似直線・曲線をScikit-learnで計算する方法についてソースコード付きでまとめました。

【例１】近似直線を回帰分析で計算

Pandasで読み込んだ信号データ(CSV)を読み込み、近似直線を回帰分析で計算します。

サンプルコード

# -*- coding: utf-8 -*-
import os
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import ntpath
from sklearn import linear_model
import joblib


# 読み込むCSVファイルのパス
csv_path = "C:/prog/python/auto/current.csv"
save_path = "C:/prog/python/"

# 空のデータフレームを作成
df = pd.DataFrame({})

# CSVファイルのロードし、データフレームへ格納
df = pd.read_csv(csv_path, encoding="UTF-8", skiprows=0)

# 電流値の列データを取り出し
Its = df.loc[:, "current"]

# 経過時間の列データを取り出し
times = df.loc[:, "time"]

# 線形回帰モデル
clf = linear_model.LinearRegression()

# 説明変数xに "x1"のデータを使用
x = np.array([times]).T

# 目的変数yに "x2"のデータを使用
y = Its.values

# 予測モデルを作成（単回帰）
clf.fit(x, y)

# パラメータ（回帰係数、切片）を抽出
[a] = clf.coef_
b = clf.intercept_

# パラメータの表示
print("回帰係数:", a)
print("切片:", b)
print("決定係数:", clf.score(x, y))

"""
回帰係数: 0.27799504539192565
切片: 147.94961709941853
決定係数: 0.3892789344510325
"""

# 学習結果の出力
joblib.dump(clf, 'C:\prog\python\scikit\clf.learn')

# 保存先パスが存在しなければ作成
if not os.path.exists(save_path):
     os.mkdir(save_path)

# グラフ化
ax = plt.axes()
plt.rcParams['font.family'] = 'Times New Roman'  # 全体のフォント
plt.rcParams['axes.linewidth'] = 1.0  # 軸の太さ

# 電流値をプロット
plt.plot(times, Its, lw=1, c="r", alpha=0.7, ms=2, label="I(t)")

# 電流値をプロット
plt.plot(times, times * a + b, lw=1, c="b", alpha=0.7, ms=2, label="Lr(t)")

# グラフの保存
plt.legend(loc="best")     # 凡例の表示（2：位置は第二象限）
plt.xlabel('Time[msec]', fontsize=12)  # x軸ラベル
plt.ylabel('Current[A]', fontsize=12)  # y軸ラベル
plt.grid() # グリッドの表示
plt.legend(loc="best")  # 凡例の表示
plt.savefig(save_path + "current.png")
plt.clf()

また、Scipyのoptimize.leastsq（最小二乗法による近似）でも同じことができます。

# -*- coding: utf-8 -*-
import os
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import ntpath
from sklearn import linear_model
import joblib
from scipy import optimize

# 1次式の近似関数
def approx_func1d(param, x, y):
    residual = y - (param[0] * x + param[1])
    return residual

# 決定係数R^2の計算
def calc_r2(x_data, y_data, y_predict):
    residuals = y_data - y_predict
    rss = np.sum(residuals**2)  # residual sum of squares = rss
    tss = np.sum((y_data - np.mean(y_data)) ** 2)  # total sum of squares = tss
    r2 = 1 - (rss / tss)
    return r2

param1 = [0, 0]

# 読み込むCSVファイルのパス
csv_path = "C:/prog/python/auto/current.csv"
save_path = "C:/prog/python/auto/"

# 空のデータフレームを作成
df = pd.DataFrame({})

# CSVファイルのロードし、データフレームへ格納
df = pd.read_csv(csv_path, encoding="UTF-8", skiprows=0)

# 電流値の列データを取り出し
Its = df.loc[:, "current"]

# 経過時間の列データを取り出し
times = df.loc[:, "time"]

# 近似直線のパラメータを計算
lq = optimize.leastsq(approx_func1d, param1, args=(times, Its))  # 係数
a = lq[0][0]
b = lq[0][1]

# 近似直線
y_predict = times * a + b

# 決定係数
r2 = calc_r2(times, Its, y_predict)

# パラメータの表示
print("回帰係数:", a)
print("切片:", b)
print("決定係数:", r2)

"""
回帰係数: 0.27799504539192565
切片: 147.94961709941853
決定係数: 0.3892789344510325
"""

# 保存先パスが存在しなければ作成
if not os.path.exists(save_path):
     os.mkdir(save_path)

# グラフ化
ax = plt.axes()
plt.rcParams['font.family'] = 'Times New Roman'  # 全体のフォント
plt.rcParams['axes.linewidth'] = 1.0  # 軸の太さ

# 電流値をプロット
plt.plot(times, Its, lw=1, c="r", alpha=0.7, ms=2, label="I(t)")

# 電流値をプロット
plt.plot(times, y_predict, lw=1, c="b", alpha=0.7, ms=2, label="Lr(t)")

# グラフの保存
plt.legend(loc="best")     # 凡例の表示（2：位置は第二象限）
plt.xlabel('Time[msec]', fontsize=12)  # x軸ラベル
plt.ylabel('Current[A]', fontsize=12)  # y軸ラベル
plt.grid() # グリッドの表示
plt.legend(loc="best")  # 凡例の表示
plt.savefig(save_path + "current.png")
plt.clf()

【例２】近似曲線を計算

Pandasで読み込んだ信号データ(CSV)を読み込み、近似曲線をScipyのoptimize.leastsq（最小二乗法による近似）で計算します。

2次の近似曲線

# -*- coding: utf-8 -*-
import os
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import ntpath
from sklearn import linear_model
import joblib
from scipy import optimize

# 近似関数
def approx_func(param, x, y):
    y = param[0]*x**2 + param[1]*x + param[2]
    return y

# 誤差関数
def residual_func(param, x, y):
    residual = y - approx_func(param, x, y)
    return residual

# 決定係数R^2の計算
def calc_r2(x_data, y_data, y_predict):
    residuals = y_data - y_predict
    rss = np.sum(residuals**2)  # residual sum of squares = rss
    tss = np.sum((y_data - np.mean(y_data)) ** 2)  # total sum of squares = tss
    r2 = 1 - (rss / tss)
    return r2


# 開始時間～終了時間
(time_start, time_end) = (0, 580)

# パラメータの初期値
param1 = [0, 0, 0]

# 読み込むCSVファイルのパス
csv_path = "C:/prog/python/auto/current.csv"

# グラフ画像の保存先パス
save_path = "C:/prog/python/auto/"

# 空のデータフレームを作成
df = pd.DataFrame({})

# CSVファイルのロードし、データフレームへ格納
df = pd.read_csv(csv_path, encoding="UTF-8", skiprows=0)

# 電流値の列データを取り出し
Its = df.loc[:, "current"]


# 経過時間の列データを取り出し
times = df.loc[:, "time"]

Its = Its[(times > time_start) & (times < time_end)]
times = times[(times > time_start) & (times < time_end)]

# 近似直線のパラメータを計算
lq = optimize.leastsq(residual_func, param1, args=(times, Its))  # 係数
(a1, a2, b) = (lq[0][0], lq[0][1], lq[0][2])

# 近似直線
y_predict = approx_func(lq[0], times, Its)

# 決定係数
r2 = calc_r2(times, Its, y_predict)

# パラメータの表示
print("回帰係数a1:", a1)
print("回帰係数a2:", a2)
print("切片:", b)
print("決定係数:", r2)

"""
回帰係数a1: -0.0011648819127047887
回帰係数a2: 1.0631165109915937
切片: 63.04535136871891
決定係数: 0.966288014400727
"""

# 保存先パスが存在しなければ作成
if not os.path.exists(save_path):
     os.mkdir(save_path)

# グラフ化
ax = plt.axes()
plt.rcParams['font.family'] = 'Times New Roman'  # 全体のフォント
plt.rcParams['axes.linewidth'] = 1.0  # 軸の太さ

# 電流値をプロット
plt.plot(times, Its, lw=1, c="r", alpha=0.7, ms=2, label="I(t)")

# 電流値をプロット
plt.plot(times, y_predict, lw=1, c="b", alpha=0.7, ms=2, label="Lr(t)")

# グラフの保存
plt.legend(loc="best")     # 凡例の表示（2：位置は第二象限）
plt.xlabel('Time[msec]', fontsize=12)  # x軸ラベル
plt.ylabel('Current[A]', fontsize=12)  # y軸ラベル
plt.grid()  # グリッドの表示
plt.legend(loc="best")  # 凡例の表示
plt.savefig(save_path + "current.png")
plt.clf()

3次の近似曲線

# -*- coding: utf-8 -*-
import os
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import ntpath
from sklearn import linear_model
import joblib
from scipy import optimize

# 近似関数
def approx_func(param, x, y):
    y = param[0] * x ** 3 + param[1] * x ** 2 + param[2] * x + param[3]
    return y

def residual_func(param, x, y):
    residual = y - approx_func(param, x, y)
    return residual

# 決定係数R^2の計算
def calc_r2(x_data, y_data, y_predict):
    residuals = y_data - y_predict
    rss = np.sum(residuals**2)  # residual sum of squares = rss
    tss = np.sum((y_data - np.mean(y_data)) ** 2)  # total sum of squares = tss
    r2 = 1 - (rss / tss)
    return r2

# 開始時間～終了時間
(time_start, time_end) = (0, 580)

# パラメータの初期値
param1 = [0, 0, 0, 0]

# 読み込むCSVファイルのパス
csv_path = "C:/prog/python/auto/current.csv"

# グラフ画像の保存先パス
save_path = "C:/prog/python/auto/"

# 空のデータフレームを作成
df = pd.DataFrame({})

# CSVファイルのロードし、データフレームへ格納
df = pd.read_csv(csv_path, encoding="UTF-8", skiprows=0)

# 電流値の列データを取り出し
Its = df.loc[:, "current"]


# 経過時間の列データを取り出し
times = df.loc[:, "time"]

Its = Its[(times > time_start) & (times < time_end)]
times = times[(times > time_start) & (times < time_end)]

# 近似直線のパラメータを計算
lq = optimize.leastsq(residual_func, param1, args=(times, Its))  # 係数
(a1, a2, a3, b) = (lq[0][0], lq[0][1], lq[0][2], lq[0][3])

# 近似直線
y_predict = approx_func(lq[0], times, Its)

# 決定係数
r2 = calc_r2(times, Its, y_predict)

# パラメータの表示
print("回帰係数a1:", a1)
print("回帰係数a2:", a2)
print("回帰係数a3:", a3)
print("切片:", b)
print("決定係数:", r2)

"""
回帰係数a1: 2.977212732750732e-06
回帰係数a2: -0.003755056988060013
回帰係数a3: 1.6641407583301977
切片: 33.97080210048596
決定係数: 0.989262775952743
"""

# 保存先パスが存在しなければ作成
if not os.path.exists(save_path):
     os.mkdir(save_path)

# グラフ化
ax = plt.axes()
plt.rcParams['font.family'] = 'Times New Roman'  # 全体のフォント
plt.rcParams['axes.linewidth'] = 1.0  # 軸の太さ

# 電流値をプロット
plt.plot(times, Its, lw=1, c="r", alpha=0.7, ms=2, label="I(t)")

# 電流値をプロット
plt.plot(times, y_predict, lw=1, c="b", alpha=0.7, ms=2, label="Lr(t)")

# グラフの保存
plt.legend(loc="best")     # 凡例の表示（2：位置は第二象限）
plt.xlabel('Time[msec]', fontsize=12)  # x軸ラベル
plt.ylabel('Current[A]', fontsize=12)  # y軸ラベル
plt.grid() # グリッドの表示
plt.legend(loc="best")  # 凡例の表示
plt.savefig(save_path + "current.png")
plt.clf()

対数関数で近似曲線

# -*- coding: utf-8 -*-
import os
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import ntpath
from sklearn import linear_model
import joblib
from scipy import optimize

# 近似関数
def approx_func(param, x, y):
    y = param[0] * np.log(x) + param[1]
    return y

# 誤差関数
def residual_func(param, x, y):
    residual = y - approx_func(param, x, y)
    return residual

# 決定係数R^2の計算
def calc_r2(x_data, y_data, y_predict):
    residuals = y_data - y_predict
    rss = np.sum(residuals**2)  # residual sum of squares = rss
    tss = np.sum((y_data - np.mean(y_data)) ** 2)  # total sum of squares = tss
    r2 = 1 - (rss / tss)
    return r2


# 開始時間～終了時間（0～500ms）
(time_start, time_end) = (0, 580)

# パラメータの初期値
param1 = [0, 0, 0]

# 読み込むCSVファイルのパス
csv_path = "C:/prog/python/auto/current.csv"

# グラフ画像の保存先パス
save_path = "C:/prog/python/auto/"

# 空のデータフレームを作成
df = pd.DataFrame({})

# CSVファイルのロードし、データフレームへ格納
df = pd.read_csv(csv_path, encoding="UTF-8", skiprows=0)

# 電流値の列データを取り出し
Its = df.loc[:, "current"]

# 経過時間の列データを取り出し
times = df.loc[:, "time"]

# 指定した時間帯のデータを抽出
Its = Its[(times > time_start) & (times < time_end)]
times = times[(times > time_start) & (times < time_end)]

# 近似直線のパラメータを計算
lq = optimize.leastsq(residual_func, param1, args=(times, Its))  # 係数
(a1, a2, b) = (lq[0][0], lq[0][1], lq[0][2])

# 近似直線
y_predict = approx_func(lq[0], times, Its)

# 決定係数
r2 = calc_r2(times, Its, y_predict)

# パラメータの表示
print("回帰係数a1:", a1)
print("回帰係数a2:", a2)
print("切片:", b)
print("決定係数:", r2)

"""
回帰係数a1: 71.40748902896294
回帰係数a2: -142.30699362830808
切片: 0.0
決定係数: 0.9611290963985782
"""

# 保存先パスが存在しなければ作成
if not os.path.exists(save_path):
     os.mkdir(save_path)

# グラフ化
ax = plt.axes()
plt.rcParams['font.family'] = 'Times New Roman'  # 全体のフォント
plt.rcParams['axes.linewidth'] = 1.0  # 軸の太さ

# 電流値をプロット
plt.plot(times, Its, lw=1, c="r", alpha=0.7, ms=2, label="I(t)")

# 電流値をプロット
plt.plot(times, y_predict, lw=1, c="b", alpha=0.7, ms=2, label="Lq(t)")

# グラフの保存
plt.legend(loc="best")     # 凡例の表示（2：位置は第二象限）
plt.xlabel('Time[msec]', fontsize=12)  # x軸ラベル
plt.ylabel('Current[A]', fontsize=12)  # y軸ラベル
plt.grid()  # グリッドの表示
plt.legend(loc="best")  # 凡例の表示
plt.savefig(save_path + "current.png")
plt.clf()