ピーク検出後その値を用いてピークフィッティングをしたいです

Question

前提

あるデータを用いて①ピーク検出②その検出した値からピークフィッティングを用いてピークを求めるのを一つにまとめたいです
色々サイトを見て進めているのでごちゃごちゃしていると思うので省略できそうな部分があればそれも教えて頂きたいです

実現したいこと

①で出たx[maxid],y[maxid]を②のx1,y1に一つずつ自動入力され値が得られるような自動解析ができるコードを作りたいです

###追記部分エラーメッセージ

RuntimeError                              Traceback (most recent call last)
C:\Users\Public\Documents\Wondershare\CreatorTemp/ipykernel_9180/2660802360.py in <module>
     88 guess_total.append(background)
     89 
---> 90 popt, pcov = curve_fit(func, x, y, p0=guess_total)
     91 
     92 

~\AppData\Local\Programs\Python\Python39\lib\site-packages\scipy\optimize\minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    792         cost = np.sum(infodict['fvec'] ** 2)
    793         if ier not in [1, 2, 3, 4]:
--> 794             raise RuntimeError("Optimal parameters not found: " + errmsg)
    795     else:
    796         # Rename maxfev (leastsq) to max_nfev (least_squares), if specified.

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 14200.

補足情報

もし必要な情報があればコメントお願いします。助けてください泣

Answer 1

①を関数にしてx1,x2を使うときに呼び出したらよいのではないですか？
理解が間違っていたらすみません。

from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import pandas as pd
import numpy as np
from scipy import signal
import os
from scipy.signal import find_peaks


x, y = np.loadtxt("./2_BLT_powder_2th-ome.txt", skiprows=29, unpack=True)


def detect_peak():
    global x
    global y
    peaks, _ = find_peaks(x, height=200)
    #signal.find_peaks(x, height=高さ, distance=距離)

    # 極大値 orderを変えることでピークの検出が変わる(ピーク検出の閾値)
    maxid = signal.argrelmax(y, order=100)
    # minid = signal.argrelmin(y, order=1) #極小値

    x_max = max(x)
    y_max = max(y)
    y_min = min(y)
    x_min = min(x)
    fig = plt.figure(figsize=(15, 5))  # 図のアスペクト比を変更（横×縦）
    Map1 = fig.add_subplot(111)
    Map1.plot(x, y)

    #plt.plot((x, y), pltsize=(6,6))
    plt.tick_params(labelsize=9)  # 目盛りの数字の大きさを変更
    plt.ylim(y_min, y_max*1.1)  # y 軸の範囲の設定,
    plt.xlim(10, 65)  # x 軸の範囲の設定
    # plt.show()

    plt.plot(x[maxid], y[maxid], 'ro', label='peak_max')
    # plt.plot(x[minid],y[minid],'bo',label='ピーク値（最小）')

    plt.legend()
    return x[maxid], y[maxid]

def func(x, *params):

    #paramsの長さでフィッティングする関数の数を判別。
    num_func = int(len(params)/3)

    #ガウス関数にそれぞれのパラメータを挿入してy_listに追加。
    y_list = []
    for i in range(num_func):
        y = np.zeros_like(x)
        param_range = list(range(3*i,3*(i+1),1))
        amp = params[int(param_range[0])]
        ctr = params[int(param_range[1])]
        wid = params[int(param_range[2])]
        y = y + amp * np.exp( -((x - ctr)/wid)**2)
        y_list.append(y)

    #y_listに入っているすべてのガウス関数を重ね合わせる。
    y_sum = np.zeros_like(x)
    for i in y_list:
        y_sum = y_sum + i

    #最後にバックグラウンドを追加。
    y_sum = y_sum + params[-1]

    return y_sum

#プロットの定義
def fit_plot(x, *params):
    num_func = int(len(params)/3)
    y_list = []
    for i in range(num_func):
        y = np.zeros_like(x)
        param_range = list(range(3*i,3*(i+1),1))
        amp = params[int(param_range[0])]
        ctr = params[int(param_range[1])]
        wid = params[int(param_range[2])]
        y = y + amp * np.exp( -((x - ctr)/wid)**2) + params[-1]
        y_list.append(y)
    return y_list

#初期値のリストを作成
#[amp,ctr,wid]

# x1, y1 = detect_peak() ##
np_x, np_y = detect_peak() ##

x2=float(input("x2"))
y2=float(input("y2"))
# guess = [] ##
guess = [[y1, x1, 1] for (x1, y1) in zip(np_x, np_y)] ##
# guess.append([y1, x1, 1]) ##
guess.append([y2, x2, 1])

#バックグラウンドの初期値
background = 70

#初期値リストの結合
guess_total = []
for i in guess:
    guess_total.extend(i)
guess_total.append(background)

popt, pcov = curve_fit(func, x, y, p0=guess_total)



fit= func(x, *popt)
plt.scatter(x, y, s=20)
plt.plot(x, fit , ls='-', c='black', lw=1)

y_list = fit_plot(x, *popt)
baseline = np.zeros_like(x) + popt[-1]
for n,i in enumerate(y_list):
    plt.fill_between(x, i, baseline, facecolor=cm.rainbow(n/len(y_list)), alpha=0.6)

print(*popt)

編集(12/01 17:35)　##が変更点です

Comments

[Nari9113]

コメントありがとうございます！①のプロット系はいらないので省略して当てはめてみましたがoperands could not be broadcast together with shapes (5328,) (23,) と出てエラーが出ました。
x[maxid]とy[maxid]が複数だからのエラーですかね、、、

[irognodyci]

完全に見落としていました、すみません。
配列扱いで修正しました。

[Nari9113]

ありがとうございます！入力してみましたところOptimal parameters not found: Number of calls to function has reached maxfev = 14200エラーが出てしまいました泣
回答して下さった部分を変更して現在のコードを質問内容の追記部分に載せさせていただきました！

[irognodyci]

どこでエラーが出ているかわかりますか？エラーメッセージと何行目かってでてくると思うのでそのまま貼り付けて質問を更新してもらえるといいかもしれません。

[Nari9113]

返信遅くなり申し訳ありません泣
エラーコード部分を質問に追加させていただきました

[irognodyci]

①でx[maxid]とy[maxid]を出力していたと思うのですが、出力されたものをそのまま全部②のほうに手動で入力したら動くんですかね？例えば、x[maxid] == [10,20,30],y[maxid] == [400,500,600]だったら、②の方で
guess.append([400,10,1])
guess.append([500,20,1])
guess.append([600,30,1])
としたら動くのかな？と。全部手動入力したら動くのであればパラメータの問題ではなく手動入力と自動入力の時とのずれになるでしょうし、それでも動かないのであればパラメータに少し手を加えないといけなくなりますし。どこに問題があるかを一つずつ整理していきたいです。

[Nari9113]

元々のコードでは仰っている②のguess.append([])を増やし手動入力は可能でした

[irognodyci]

となると、うまく入ってないんですかね。
手動入力後にprint(guess)した場合と、自動入力の方で入力した後print(guess)した場合とで違いはありますか？

[Nari9113]

自動入力で
[[131.0, 10.755, 1], [63.0, 12.205, 1], [58.0, 13.395, 1], [262.0, 16.195, 1], [59.0, 18.395, 1], [194.0, 21.665, 1], [272.0, 23.355, 1], [72.0, 24.635, 1], [113.0, 26.925, 1], [1286.0, 30.055, 1], [440.0, 33.085, 1], [43.0, 34.725, 1], [88.0, 37.025, 1], [251.0, 38.365, 1], [240.0, 39.775, 1], [77.0, 42.855, 1], [287.0, 47.445, 1], [256.0, 51.445, 1], [86.0, 53.505, 1], [39.0, 55.445, 1], [329.0, 57.245, 1], [71.0, 58.555, 1], [134.0, 62.535, 1]]
と出たので違いはなさそうです

[irognodyci]

元ファイルがないのでこちらで検証のしようがないのが痛いところです。原因がすぐわからずいろいろさせてしまってすみません、、
手動入力してprintしたものをコピーして、自動入力の方のプログラムで、guess = [[y1, x1, 1] for (x1, y1) in zip(np_x, np_y)]をguess = printしたもの(コピーしたものを貼り付ける)
にしても動作しますか？

[Nari9113]

実験結果のファイルなので渡したいぐらいです、、いえいえ！進んでる気がしてむしろ喜んでいます笑
ちゃんと理解して入力できているか分かりませんがguess = [[1286, 30.055, 1],[440, 33, 1]] のように入れ込めば表示されました

[Nari9113]

計算が多すぎる問題ですかね、、(´；ω；｀)

[Nari9113]

このデータ、全部入れてみようとしたらできず、後ろ4組分消したところ反映されました
ピークのノイズ分も計算して後から消そうとしていたので計算量が多くなったのかと泣
最初ノイズ部分を消せなかったのでノイズ部分を消せたら全部動くかもしれません泣

[irognodyci]

手動入力してたときは少ないデータでやっていたってことですか？
ノイズまで拾っていると思うと全部手動入力は厳しそうですが(-_-;)

[Nari9113]

まずプログラムが動くかを試していたので少ないデータで行っていました、、
ある程度のノイズ除去まではしていたので、、まだ多かったみたいです泣

[Nari9113]

guessの入れ込む前にノイズと思われるyがバックグラウンドの二倍ほどであればその組を省くということは可能ですか？

[irognodyci]

background = 70
guess = [[y1, x1, 1] for (x1, y1) in zip(np_x, np_y) if background * 2 > y1]
でできると思います。backgroundを先に定義したのと、リスト内包表記の後ろの方にif文付けただけです！

[irognodyci]

そもそもノイズをピークとして拾ってる場合ってピークフィッティングできるんでしょうか？

[Nari9113]

ありがとうございます！
ノイズを除いたデータ(元のデータとは別)で行っていたので可能でした
最初からまとめて行うと訳が分からなくなりそうだったので、、

[Nari9113]

バックグラウンド2倍以下除去の関数を打ち込んだところエラーコードの最後がmaxfev = 6400に変わりました

[irognodyci]

popt, pcov = curve_fit(func, x, y, p0=guess_total, maxfev=14200)
maxfevを変えると推測回数が変わるので、すごく大きくしたらいける、、のか？
この辺は触ったことないので僕もわかってなくて

[Nari9113]

50000も駄目で100000にしたらいけました！！！！！！(いや関数呼び出し多すぎ、、)
後は出た値をリスト化するだけになりました！本当にありがとうございます！！！

[irognodyci]

それでいいのか感は少しありますが、とりあえず動いて良かったです！