Go to Qiita Advent Calendar 2024 Top
LoginSignup
Tetsu27
@Tetsu27

Are you sure you want to delete the question?

Leaving a resolved question undeleted may help others!

定式化を変えてフィッテングしたい

Q&A

Closed

Pythonscipyfitting

解決したいこと

横軸によってスケールの異なるデータ、例えば下図、をフィッテングする際に次の定式化でフィッテングしたいです。
$y_{i}$を観測量、$\hat{y_{i}}$をフィッティングカーブ上の$x_i$での量、$\vec{a}$をフィッティングカーブのパラメータを並べたベクトルとするときに

\text{min} \left( \sum_{i=1}^{m} \frac{y_i - \hat{y_i}(x_i; \vec{a})}{y_i} \right)

でフィッティングをしたいです。
image.png

つまり、観測量のスケールが大きなところは重要視させず、信号の観測量の小さいところを重要視させフィッティングさせるということがしたいです。

これはscipyのcurve_fitで実現できるものなのでしょうか。色々と考えたのですがフィッティングアルゴリズム自体の変更なので、curve_fitでは実装は難しいと考えておりました。
しかし、何かアイデアや実装方法ご存じの方がいましたらお教えいただきたいです。
よろしくお願いいたします。。

0

1Answer

PondVillege
@PondVillege

色々と考えたのですがフィッティングアルゴリズム自体の変更なので、curve_fitでは実装は難しいと考えておりました。

そうですね，質問者の用意した数式を用いて

観測量のスケールが大きなところは重要視させず、信号の観測量の小さいところを重要視させフィッティングさせるということがしたい

というような，モデルに制約を課すことはよくある話ですが，curve_fitでは既存アルゴリズムしか選択できず，質問者の使いたい制約を課すためには，ライブラリ自体を直接いじるほかありません．

自分の知っている範囲だと，例えば機械学習系ライブラリを使った関数フィッティングを推奨しておきます．

$\text{sinc}(x)=\frac{\sin(ax+b)}{ax+b}$の$a,b$を求めるべく，次のように書くことができました．

コード 
import tensorflow as tf
from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau
import numpy as np

# ゼロ除算を避けるための小さな値
eps = 1e-10
# sinc関数のa, b
a, b = 2, 3

def f(x, a, b):
    x = a * x + b
    return tf.sin(x) / (x + eps)
    
class ParametricSincModel(tf.keras.Model):
    def __init__(self, lmbd=0.001):
        super(ParametricSincModel, self).__init__()
        # `a`を1.0、`b`を0.0で初期化
        self.a = tf.Variable(initial_value=1.0, dtype=tf.float32, trainable=True)
        self.b = tf.Variable(initial_value=0.0, dtype=tf.float32, trainable=True)
        self.lmbd = lmbd # lambda: custom_loss制約の強さ

    def call(self, inputs):
        # フィットしたい関数fを利用
        return f(inputs, self.a, self.b)

    def custom_loss(self, y_true, y_pred):
        # 損失関数に(y_true - y_pred) / (y_true + eps)を追加して、y_trueが小さい時に大きな重みがかかるようにする
        return tf.reduce_mean(tf.abs((y_true - y_pred) / (y_true + eps)))

    def train_step(self, data):
        x, y_true = data
        with tf.GradientTape() as tape:
            y_pred = self.call(x)
            # 既存の損失（この例ではMSE）を計算
            loss = self.compiled_loss(y_true, y_pred, regularization_losses=self.losses)
            # カスタム損失を計算
            custom_loss = self.custom_loss(y_true, y_pred)
            # 既存の損失にカスタム損失を加算
            total_loss = loss + self.lmbd * custom_loss
        gradients = tape.gradient(total_loss, self.trainable_variables)
        self.optimizer.apply_gradients(zip(gradients, self.trainable_variables))
        # 更新された損失を返す
        return {'loss': loss, 'custom_loss': custom_loss, 'total_loss': total_loss}

# sinc関数からデータを生成
def generate_data(num_samples=100000):
    x = tf.linspace(-10.0, 10.0, num_samples)
    y = f(x, a, b) # 上で定義したa, bを使う．観測値が描いている関数はこの関数
    return x, y

# データを生成
x_train, y_train = generate_data()
# 入力データの形状を調整（モデルに合わせて）
x_train = tf.reshape(x_train, (-1, 1))

# モデルのインスタンスを作成
model = ParametricSincModel()

# コンパイル：オプティマイザと損失関数を設定
model.compile(optimizer='adam', loss='mse')

# EarlyStoppingの設定：指定されたpatience回のエポック数以上、
# val_lossが改善されなければトレーニングを停止する
early_stopping = EarlyStopping(
    monitor='total_loss',
    patience=10,
    restore_best_weights=True
)

# ReduceLROnPlateauの設定：指定されたpatience回のエポック数以上、
# val_lossが改善されなければ学習率を減少させる
reduce_lr = ReduceLROnPlateau(
    monitor='total_loss',
    factor=0.1,
    patience=5,
)

# トレーニング：モデルをデータにフィットさせる
model.fit(
    x_train,
    y_train,
    epochs=1000,
    callbacks=[early_stopping, reduce_lr], 
    #verbose=0
)

# 学習したパラメータ a と b の値を取得
a_value, b_value = model.a.numpy(), model.b.numpy()
print(f"Learned parameters: a = {a_value}, b = {b_value}")
print(f"Absolute difference: a = {np.abs(a_value - a)}, b = {np.abs(b_value - b)}")
標準出力 
Epoch 1/1000
2048/2048 [==============================] - 1s 459us/step - loss: 0.0902 - custom_loss: 10.2247 - total_loss: 1.1491 - lr: 0.0010
Epoch 2/1000
2048/2048 [==============================] - 1s 438us/step - loss: 0.1754 - custom_loss: 10.5889 - total_loss: 1.1869 - lr: 0.0010
Epoch 3/1000
2048/2048 [==============================] - 1s 449us/step - loss: 0.0177 - custom_loss: 7.4664 - total_loss: 0.8766 - lr: 0.0010
Epoch 4/1000
2048/2048 [==============================] - 1s 440us/step - loss: 0.1671 - custom_loss: 9.4329 - total_loss: 1.0701 - lr: 0.0010
Epoch 5/1000
2048/2048 [==============================] - 1s 443us/step - loss: 0.0952 - custom_loss: 7.5609 - total_loss: 0.8857 - lr: 0.0010
Epoch 6/1000
2048/2048 [==============================] - 1s 477us/step - loss: 0.0873 - custom_loss: 10.2695 - total_loss: 1.1589 - lr: 0.0010
Epoch 7/1000
2048/2048 [==============================] - 1s 459us/step - loss: 0.1193 - custom_loss: 9.0739 - total_loss: 1.0395 - lr: 0.0010
Epoch 8/1000
2048/2048 [==============================] - 1s 444us/step - loss: 0.1108 - custom_loss: 10.5612 - total_loss: 1.1862 - lr: 0.0010
Epoch 9/1000
2048/2048 [==============================] - 1s 448us/step - loss: 0.1211 - custom_loss: 10.4570 - total_loss: 1.1757 - lr: 1.0000e-04
Epoch 10/1000
2048/2048 [==============================] - 1s 443us/step - loss: 0.1399 - custom_loss: 10.5988 - total_loss: 1.1900 - lr: 1.0000e-04
Epoch 11/1000
2048/2048 [==============================] - 1s 444us/step - loss: 0.1457 - custom_loss: 10.0044 - total_loss: 1.1310 - lr: 1.0000e-04
Epoch 12/1000
2048/2048 [==============================] - 1s 446us/step - loss: 0.1175 - custom_loss: 13.9249 - total_loss: 1.5210 - lr: 1.0000e-04
Epoch 13/1000
2048/2048 [==============================] - 1s 452us/step - loss: 0.1228 - custom_loss: 19.6584 - total_loss: 2.0963 - lr: 1.0000e-04
Learned parameters: a = 1.0008928775787354, b = -0.06013570353388786
Absolute difference: a = 0.9991071224212646, b = 3.060135703533888

このようになり，ある程度の成功を収めているように見えます．
しかし，curve_fitで使われるような二乗誤差の最小化に加えて，質問者が使いたい制約関数custom_lossを付与してしまうと，精度が下がることがわかりました．実際，コード中でtotal_loss = loss + self.lmbd * custom_lossと書いてますが，このself.lmbd=0.001を1まで増やした場合，実際の$a,b$とは遠い値になってしまいました．self.lmbd=0のときに，$a, b = 2, 3$ちょうどになります．つまり，custom_lossはない方がいいということです．

次のself.lmbd=0.1のときのグラフを見てください

output.png

ここでは，$y=0$付近の制約が強すぎて，ここさえ一致していれば良い．みたいな局所最適解に陥った例です．この問題はcurve_fitでも同様で初期値依存の解を出すため，これはcustom_lossを用意してかつ初期値$a,b=1,0$を設定したことが原因と言えます．実際に，初期値からあまり動かず，最終的な解は$a,b=1.001,-0.06$となっています．$a,b=2,3$まで動かずして局所最適解に陥ってますね．

具体的には，標準出力を見てわかるようにcustom_lossの値が，MSEで得られたlossに比べて非常に大きいのが原因です．これにつられて$a,b$の最適化を進めることになってしまい，結果として局所最適解を出しました．

しかしこれはあくまで$\text{sinc}(x)=\frac{\sin(ax+b)}{ax+b}$の$a,b$を求めるべくして起こった事象ですので，質問者のデータの範囲でやってみながら検証するのが良いと思います．

1Like

Comments

  1. @Tetsu27

    Questioner

    返信が遅くなり申し訳ありません。通知メールに気が付かず返信が遅くなってしまいました。
    大変丁寧なご回答誠にありがとうございます。
    curve_fitではライブラリをいじるほか、対策はない旨理解しました。
    機械学習系ライブラリを使った関数フィッティングをコード付きで、お教えいただき本当にありがとうございます。局所最適解に陥る危険性がある旨、理解致しました。こちらを自分のデータで挑戦してみます。
    この度は大変丁寧にお教えいただき、本当にありがとうございます。大変助かりました。

  2. @PondVillege

    ちなみに観測データみたくノイズを含ませてフィッティングさせたところ，self.lmbd=0.1でも局所最適解に陥らず大域的最適解を得ることができました．独自の制約を使わなくても局所最適解に陥る例を見てきたので，質問者がうまくいくことを祈っております．

Your answer might help someone💌

LoginSign Up