pythonで曲線の各点の傾きを表すグラフを作成したいです

Question

前提・実現したいこと

pythonで曲線の各点の傾きを表すグラフを作成したいです。
言い換えるならば、x-y の配列・グラフから、x-y’の配列・グラフを求める方法が知りたいということです。
自分なりに配列の微分や勾配の求め方を調べ、numpyのnp.gradientという機能を用いて作成したのですが、x-yの曲線から想像できるような滑らかな曲線とはなりませんでした。
他によりよい方法をご存じの方がいらっしゃいましたらご教授いただけると幸いです。
（今回解析しているのは磁化特性曲線で、横軸の磁界に対する縦軸の磁束の変化率（透磁率）を求めようとしています。）

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt

x = H[209:416] #磁界の値をｘに代入
y = B[209:416] #磁束の値をｙに代入

plt.figure()
plt.plot(x,y)
plt.grid()

plt.show()

各点の傾きを求めたい曲線(画像1)

dx=np.gradient(x) #xの微小分
dy=np.gradient(y) #yの微小分
myu=dy/dx　　　　 #各点の微分した値
plt.figure()
plt.plot(x,myu)
plt.grid()

plt.show()

試作したグラフ(画像2)

発生している問題・エラーメッセージ

(画像1)で線形となっている部分が(画像2)のようにギザギザになってしまっています。より適した微分の方法などをご存知でしたら回答よろしくおねがいします。

Answer 1

【追記】
xの間隔（座標）を考慮すると質問者さんと同じようなグラフになりました。

dy = np.gradient(y, x)
plt.plot(x, dy)
plt.show()

---

【xの間隔が1の場合】

dy = np.gradient(y)
plt.plot(x, dy)
plt.show()

【NumPy】配列の勾配 numpy.gradient

Comments

[KNKR]

ご回答ありがとうございました！np.gradientは差分を見ているのではなく、微分の定義そのもののようなコマンドだということですね。助かりました！

[meg_]

すみません。上記もxの間隔がデフォルト値の1の場合でした。

[meg_]

xの差分を考慮しようと思ったのですが、良く見るとxの並びが前後している部分がいくつかありますね。
例：397.75174812, 397.69494205　←　xの後ろから3番目と2番目
そのため計算できませんでした。

[KNKR]

ちょうど質問しようとしていたところでした(笑）
やはり間隔がバラバラだと使えないですかね…。
もし、xの並びが前後しないような綺麗なデータが取れた時はどうするのがいいでしょうか。お時間ありましたらコメントお願いいたします。

[KNKR]

追記ありがとうございました！また試行錯誤してみます！

[meg_]

データが関数で作成されたものではなく実測値の場合はこのようなプロットになりませんか？(ばらつきの範囲とは考えられませんか？)

[KNKR]

考えられますね。やはり実測値の微小分を用いた微分では仕方ないですよね。
測定方法などを見直してみます。対応本当にありがとうございました。

Answer 2

Savitzky–Golay法

# 平滑化点数5
dy5 = np.convolve(y, np.arange(2, -3, -1), 'varid') / 10
myu5 = dy5 / dx[2:-2]

# 平滑化点数7
dy7 = np.convolve(y, np.arange(3, -4, -1), 'varid') / 28
myu7 = dy7 / dx[3:-3]

# 平滑化点数9
dy9 = np.convolve(y, np.arange(4, -5, -1), 'varid') / 60
myu9 = dy9 / dx[4:-4]

plt.plot(x, myu, x[2:-2], myu5, x[3:-3], myu7, x[4:-4], myu9)
plt.grid()
plt.show()

点数が多いほど平滑化が強くなります

【訂正】
これはxが等間隔じゃないと使えない手法です
データ見たら等間隔ではないので、使えませんね
失礼しました

【追記】
Savitzky-Golay smoothing filter for not equally spaced data
に、不等間隔(non-uniform)でも計算できるSavitzky–Golay法のPythonコードが載ってたので、それを微分(differential)にも使えるように修正しました
よろしければお試しください

import numpy as np
import math

def non_uniform_savgol_der(x, y, window, polynom, der):
    """
    Applies a Savitzky-Golay filter to y with non-uniform spacing
    as defined in x

    This is based on https://dsp.stackexchange.com/questions/1676/savitzky-golay-smoothing-filter-for-not-equally-spaced-data
    The borders are interpolated like scipy.signal.savgol_filter would do

    Parameters
    ----------
    x : array_like
        List of floats representing the x values of the data
    y : array_like
        List of floats representing the y values. Must have same length
        as x
    window : int (odd)
        Window length of datapoints. Must be odd and smaller than x
    polynom : int
        The order of polynom used. Must be smaller than the window size
    der : int
        The order of derivative. Must be positive and smaller than polynom order

    Returns
    -------
    np.array of float
        The smoothed y values
    """
    if len(x) != len(y):
        raise ValueError('"x" and "y" must be of the same size')

    if len(x) < window:
        raise ValueError('The data size must be larger than the window size')

    if type(window) is not int:
        raise TypeError('"window" must be an integer')

    if window % 2 == 0:
        raise ValueError('The "window" must be an odd integer')

    if type(polynom) is not int:
        raise TypeError('"polynom" must be an integer')

    if polynom >= window:
        raise ValueError('"polynom" must be less than "window"')

    if polynom < 0:
        raise ValueError('"polynom" must be larger than 0')

    if type(der) is not int:
        raise TypeError('"der" must be an integer')

    if der < 0:
        raise TypeError('"der" must be an positive integer')

    if der > polynom:
        raise TypeError('"der" must be equal or less than "polynom"')

    half_window = window // 2
    polynom += 1

    # Initialize variables
    A = np.empty((window, polynom))     # Matrix
    tA = np.empty((polynom, window))    # Transposed matrix
    t = np.empty(window)                # Local x variables
    y_smoothed = np.full(len(y), np.nan)

    # Start smoothing
    for i in range(half_window, len(x) - half_window, 1):
        # Center a window of x values on x[i]
        for j in range(0, window, 1):
            t[j] = x[i + j - half_window] - x[i]

        # Create the initial matrix A and its transposed form tA
        for j in range(0, window, 1):
            r = 1.0
            for k in range(0, polynom, 1):
                A[j, k] = r
                tA[k, j] = r
                r *= t[j]

        # Multiply the two matrices
        tAA = np.matmul(tA, A)

        # Invert the product of the matrices
        tAA = np.linalg.inv(tAA)

        # Calculate the pseudoinverse of the design matrix
        coeffs = np.matmul(tAA, tA)

        # Calculate c0 which is also the y value for y[i]
        y_smoothed[i] = 0
        for j in range(0, window, 1):
            #y_smoothed[i] += coeffs[0, j] * y[i + j - half_window]
            y_smoothed[i] += coeffs[der, j] * y[i + j - half_window] * math.factorial(der)

        # If at the end or beginning, store all coefficients for the polynom
        if i == half_window:
            first_coeffs = np.zeros(polynom)
            for j in range(0, window, 1):
                for k in range(polynom):
                    first_coeffs[k] += coeffs[k, j] * y[j]
        elif i == len(x) - half_window - 1:
            last_coeffs = np.zeros(polynom)
            for j in range(0, window, 1):
                for k in range(polynom):
                    last_coeffs[k] += coeffs[k, j] * y[len(y) - window + j]

    # Interpolate the result at the left border
    for i in range(0, half_window, 1):
        y_smoothed[i] = 0
        x_i = 1
        #for j in range(0, polynom, 1):
        for j in range(der, polynom, 1):
            y_smoothed[i] += first_coeffs[j] * x_i * math.factorial(j)
            x_i *= x[i] - x[half_window]

    # Interpolate the result at the right border
    for i in range(len(x) - half_window, len(x), 1):
        y_smoothed[i] = 0
        x_i = 1
        #for j in range(0, polynom, 1):
        for j in range(der, polynom, 1):
            y_smoothed[i] += last_coeffs[j] * x_i * math.factorial(j)
            x_i *= x[i] - x[-half_window - 1]

    return y_smoothed

if __name__ == '__main__':
    import numpy as np
    import matplotlib.pyplot as plt

    # テスト用ダミーデータ
    x = np.arange(0, 10, 0.1)
    y = 2 * x * x + 3 * x + 4 + np.random.randn(100) / 10
    y11 = non_uniform_savgol_der(x, y, 11, 2, 0)
    plt.plot(x, y, x, y11)
    plt.grid()
    plt.show()

    # 微分
    dy = np.gradient(y, x)
    dy11 = non_uniform_savgol_der(x, y, 11, 2, 1)
    plt.plot(x, dy, x, dy11)
    plt.grid()
    plt.show()

    # 二階微分
    ddy = np.gradient(dy, x)
    ddy11 = non_uniform_savgol_der(x, y, 11, 2, 2)
    plt.plot(x, ddy, x, ddy11)
    plt.grid()
    plt.show()

xが順番に並んでないのは、sortすれば直せます

# テスト用ダミーデータ
x = np.array([1, 2, 3, 5, 4, 6])
y = x * 2
print(x)
print(y)

# ソート (xの順番でyも)
tmp = zip(x, y)
tmp2 = sorted(tmp)
xx, yy = zip(*tmp2)
xx = np.array(xx)
yy = np.array(yy)

# 比較
print(x)
print(xx)
print(y)
print(yy)

Comments

[KNKR]

なめらかな曲線が描けました。対応本当にありがとうございました。

[jbpb0]

関数の三つ目の引数(window)を大きくするほど、平滑化が強くなり細かいデコボコが減りますが、全体の形状(左右の急激な立ち上がり・立ち下がり等)が徐々に崩れてきます
そのバランスで数値を決めますので、最適な数値はデータや目的によってケースバイケースです

[jbpb0]

二階微分の処理にバグがあったので、関数non_uniform_savgol_derのコードを差し替えました
なお、(二階ではない単なる)微分には影響ありません