数値計算の結果から測定値にフィットするパラメータを探索したい

Question

前提・実現したいこと

研究用途として、数値計算のプログラムを作成しています。
解析的には解くことができない偏微分方程式の数値計算の結果として、下の表のように時間ステップt=0, 1, 2,..., iごとの計算結果yのリストが格納されたresult.csvファイルが出力されています。

t＼x	0	1	...	20
0	y[0,0]	y[0,1]	...	y[0,20]
1	y[1,0]	y[1,1]	...	y[1,20]
2	y[2,0]	y[2,1]	...	y[2,20]
...
i	y[i,0]	y[i,1]	...	y[i,20]

一方、実データdata.csvは下記のような形式になっています。

x	0.00	0.74	1.08	...	19.39
y	37	24	59	...	43

この結果から、実データを最も適切に説明する時間tを探索する方法を実装しようとしています。

発生している問題

例えば以下のようにして、時間ステップごとの測定値-数値計算の結果の残差二乗和を計算し、その最小値をmin()で探せば良いと考えました。

import numpy as np
import pandas as pd

df_result = pd.read_csv("result.csv", index_col=0).values
df_data = pd.read_csv("data.csv", index_col=0).iloc[0].values
list_rsds = []
for i in range(df_result.shape[0]):
    rsd_sum = 0
    for j in range(len(df_data)):
        rsd_sq = (df_data[j] - df_result[i, j])**2
        rsd_sum += rsd_sq
    list_rsds.append(rsd_sum)
print(list_rsds.index(min(list_rsds)))

しかし、実データのx座標の値が等間隔ではないため、数値計算の結果と実データのx座標が異なり（定数倍ですらない）、直接残差を計算することができません。

試したこと

「パラメータフィッティング python」等で検索しましたが、フィッティングする関数が予め与えられている場合がヒットしました。
補間やカーブフィッティングなどの最適化
しかし、数値計算の結果のような離散値をどのように扱うべきか分かりませんでした。

解決案として、

1. 測定値を線形補間して2点間の直線の式を求め、数値計算のx座標の値を代入->残差を計算する
2. 数値計算の値を何らかの多変数関数として表現し、その関数を用いてパラメータフィッティングを行う

~~
といったことを考えているのですが、このような方針がそもそも適切なのか、どのように実装するのが良いのか方針が立っておりません。~~
1で示した方法についてコードを書いてみましたが、もう少しうまいやり方があるような気がするのでアドバイスを頂けると幸いです。

# 直線の式y = a * x + bのa, bを計算してparamsとして出力
def calc_param(result):
    dx = 1 / len(result)
    result_distance = np.linspace(0,1,len(result),endpoint=True) # x座標を0-1に規格化している
    params = np.zeros([len(result)-1, 2])
    for j in range(len(result)-1):
        a = (result[j+1] - result[j]) / dx
        b = result[j] - j * dx
        params[j] = [a,b]
    return params, result_x

def main():
    df_result = pd.read_csv("result.csv", index_col=0).values

    data_value = [37, 24, 59,..., 43]
    data_x = [0, 0.74, 1.08,..., 19.39]

    # x座標の規格化
    data_x = data_x / max(data_x)

    # 計算
    list_reses = []
    for i in range(len(df_result)):
        result = df_result[i]
        params, result_x = calc_param(result)
        reses = 0
        #print(len(result_x))
        # 数値計算のx座標を順番に入力
        for m in range(len(result_x)-1):
            # 分析値のx座標を順番に入力
            for n in range(len(data_x)):
                # もし範囲内にあれば内挿して残差二乗を計算
                # 末尾は規格化して1なのでresult_x[m+1]==data_x[n]
                if ( result_x[m] <= data_x[n] < result_x[m+1] ) \
                   or (result_x[m+1] == data_x[n]):
                    f = params[m,0] * data_x[n] + params[m,1]
                    res = (f - data_value[n]) ** 2
                    reses += res
                else:
                    pass
        list_reses.append(reses)
    print(list_reses.index(min(list_reses)))
if __name__ == '__main__':
    main()

補足情報（FW/ツールのバージョンなど）

• OS: Windows 10 Home
• Python 3.8.5

Comments

[hayataka2049]

研究ですと、数値計算側で測定データのxに合わせてデータを吐くか、数値計算のxに合わせて測定するかが一番無難では？　論文の書きやすさで言えば。

[退会済みユーザー]

天然物を扱っているため、測定点を等間隔に取るのが難しく、このような形になっています。数値計算の際の⊿xを測定点に合わせて可変にすることも検討いたしましたが、計算量が多くなってしまうためフィッティングの方法を工夫する方向でできないか考えています。

Answer 1

とりあえず，「差」をy値の差だけで見て良い話なのだとすれば，
実データのx値に

偏微分方程式の数値計算の結果

のx値を合わせればよいだけに思えます．
そうしない理由が全くわかりません．

---

「何故かxが実データと合わされていない数値計算の結果」　から　
「実データの各x値におけるy値」をどうにかして用意してやる，のだとして，

そのための方法が
スプライン補間なのか多項式当てはめ等になるのか，その他なのか…？
どんな方法が「適切」のかは扱っている関数の意味や形状に依存するでしょうけども，
いずれにせよ，
そのようにして適当に値を補間(?)してやること自体の正当性に関して十分な説得力が必要に思います．

Comments

[fana]

xが1刻みというのが，十分に「細かい」データなのであって，
x=1 と x=0.74 や x=1.08 との差というのが十分に「微小」と言えるスケール感（？）なのかな…？

[退会済みユーザー]

回答ありがとうございます。コメントの認識の通りです。
前者のアプローチを取れないのは実データのx値が等間隔には引けないためです。
後者に関しては線形補完が適切だろうと考えています。

[fana]

これ以上細かく数値計算すること　と　補間すること　との間にはもはや意味のある差はない，と言える状態なら，補間しても良いのかも．

[fana]

（ここでは「数値計算」自体の精度には問題が無いのだとして）

もう少し刻み幅を小さくして（1/2とか1/4とか？）やってみたやつと，現状の1刻みの値から補間した値との差くらいは見ておいた方がよいかもしれません．念のため．

[退会済みユーザー]

アドバイスありがとうございます。詳細はここに書けないのですが、少なく見積もって1σ=0.5以上は誤差がつくので測定値も計算値も十分に細いだろうとラボ内では考えています。より細かくした場合については時間ステップとの兼ね合いもありますが検討したいと思います。

Answer 2

scipyなどを使えば離散値のカーブフィッティング自体はなんとでも出来そうですが、フィッティングすること自体の意味については、研究内容によると思うので、同じ研究室の方などと相談される方がいいと思います。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
https://qiita.com/maskot1977/items/e4f5f71200180865986e

Comments

[退会済みユーザー]

回答ありがとうございます。（質問内容がうまく伝わっていないか私の知識不足でScipyを誤解しているかもしれないので確認したいのですが)「試したこと」の項で上げた2つの解決策のうち、前者（測定値を線形補間）はあまり現実的ではないということでしょうか？ドキュメントを読む限りでは、何らかの関数を用意する必要があるように思われるので。数値計算の内容は解析的には解けない偏微分方程式になっていて、関数の形で表すことはできません。※わかりにくい文章になっていたため追記しました。すみません

[hide5stm]

フィッティング自体は、線形でも非線形でも、scipyなどを駆使すればいくらでも出来ると思います。
データが何関数でフィッティングされるべきかは、研究内容がわかる人しかわからないですし、
プログラミングの領域からも外れていきます。

[退会済みユーザー]

ありがとうございます。Scipyは使ったことがないので、ドキュメントをもう少し読んで考えてみます。もしフィッティングをする場合は、どのような関数を使うべきかに関してはこちらでしっかりと判断したいと思います。
※整理がついておらず、質問がピンぼけになっていました。線形補間による一応の案を追記いたしました(当初の質問内容からかなり外れてしまったので、お忙しければスルーしていただいて結構です)

Answer 3

追記したコードで「1. 線形補間した上で残差を計算」は達成できるので自己解決とさせていただきます