Pythonでのガウシアンフィッティングでエラーが出る

Question

Pythonでガウシアンフィッティングしようとすると以下のようなエラーが出てしまいます。

 File "/Users/pyworks/gaus.py", line 54, in <module>
    paramater_optimal, covariance = curve_fit(fitting, channels, counts, initial_parameter)
  File "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/scipy/optimize/minpack.py", line 789, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
  File "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/scipy/optimize/minpack.py", line 410, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/scipy/optimize/minpack.py", line 24, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/scipy/optimize/minpack.py", line 485, in func_wrapped
    return func(xdata, *params) - ydata
ValueError: operands could not be broadcast together with shapes (1500,) (1500,3) 

実行したコードは以下です。
https://lastline.hatenablog.com/entry/2018/04/12/112502
を参考にしています。このサイトに載ってるコードのフィッティングする元データと初期値のみいじっております。

#モジュールのインポート
import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import cv2

#フィッティングしたいデータの生成
img = cv2.imread('convert.png') #1500×360の画像読み込み
img_split = np.split(img, 36)
#元々の配列(1500×360)の時、縦を36分割、1500×10の配列を36個生成

#フィッティングしたいデータ
channels=list(range(1500))
counts=np.sum(img_split[0],axis=0)


#フィッティング関数の定義
def quadratic(x, a, b, c):
    return a + b * x + c * x**2

#ガウシアンはより簡素に記述できるが、初期値を類推しやすい式にした。
def gauss(x, a, sigma, mean):
    return a / np.sqrt(2.0*np.pi) / sigma * np.exp(-((x-mean)/sigma)**2/2)

#今回はあまりにも冗長になるため、*paramaters と指定した。
def fitting(x, *paramaters):
    a0, b0, c0, a1, sigma1, mean1, a2, sigma2, mean2, a3, sigma3, mean3, = paramaters
    return quadratic(x, a0, b0, c0) + \
           gauss(x, a1, sigma1, mean1) + \
           gauss(x, a2, sigma2, mean2) + \
           gauss(x, a3, sigma3, mean3)

#初期値の設定。
#ガウシアンの場合、meanを適切に設定しないと発散しやすい。
#線形近似であっても適切に設定することが望ましい。
#関数の引数が複数ある場合は、タプルで渡すこと。
initial_quadratic = 0.1, 0.1, 0.1 #a0, b0, c0
initial_gauss1 = 2500, 464, 170. #a1, sigma1, mean1
initial_gauss2 = 2400, 1000, 100. #a2, sigma2, mean2
initial_gauss3 = 2000, 1300, 100. #a3, sigma3, mean3

#初期値が多いので、個別に指定して接続した。
initial_parameter = initial_quadratic + initial_gauss1 + initial_gauss2 + initial_gauss3

#フィッティングの実行
#paramater_optimalが最適化されたパラメータで、covarianceは共分散。
paramater_optimal, covariance = curve_fit(fitting, channels, counts, initial_parameter)

#matplotlibの設定
plt.rc('font', family='Arial', size=14) #フォントの設定
plt.rc('xtick.major', width=1, size=6) #x軸の主目盛りの設定
plt.rc('xtick', direction='in', top=True) #x軸目盛りの向き、上側に表示するか
plt.rc('ytick.major', width=1, size=6) #y軸の主目盛りの設定
plt.rc('ytick', direction='in', right=True) #y軸目盛りの向き、右側に表示するか
plt.rc('axes', linewidth=1.5) #枠線の太さ
plt.rc('lines', linewidth=1.0) #プロットの太さ

#今回は、Axes が一つしかないので fig, ax1 = plt.subplots() の方が簡素。
#慣例として ax1 を用いるが、df 同様に意味のある変数名にすべき。
fig = plt.figure(figsize=(8, 6), dpi=80)
gs = gridspec.GridSpec(1, 1)
ax1 = fig.add_subplot(gs[0])

#データをグラフに表示する。
ax1.plot(channels, counts, ls='-', c='g', label='row data')
ax1.set_xlabel('Channel', fontsize=16)
ax1.set_ylabel('Count', fontsize=16)
ax1.yaxis.set_data_interval(0, 100) #set_ylim を使う方が一般的

#フィッティングした曲線の表示
#np.corrcoef で相関係数を求める。
#結果が2×2の行列で返るので、1行目2列目もしくは2行目1列目の値を取得。
ax1.plot(channels, fitting(channels, *paramater_optimal),
         lw=1.5, ls='--', c='k', label='fitting curve\n' + \
         f'$R^2$ = {np.corrcoef(counts, fitting(channels, *paramater_optimal))[0][1]}')

#二次関数のパラメーターを利用したいので、スライスで取り出す。
a0, b0, c0 = paramater_optimal[0:3]

ax1.plot(channels, quadratic(channels, a0, b0, c0), c='b',
         label=r'$y = a_0 + b_0x + c_0x^2$' +'\n'+ f'$a_0$ = {a0}\nb$_0$ = {b0}\nc$_0$ = {c0}')

#各meanをグラフ上に表示する。繰り返しのため for 文を利用。
for i in range(3,12,3):
#この関数内の、a0, b0, c0 はグローバル変数である。
    def quadratic_gauss(x, *paramaters):
        a1, sigma1, mean1, = paramaters
        return quadratic(x, a0, b0, c0) + \
               gauss(x, a1, sigma1, mean1)

#annotate で Axes 内にテキストや矢印を表示できる。
    mean = paramater_optimal[i+2]
    ax1.annotate(int(mean), fontsize = 14,
                 xy=(mean, quadratic_gauss(mean, *paramater_optimal[i:i+3])+3),
                 xytext=(mean, quadratic_gauss(mean, *paramater_optimal[i:i+3])+15),
                 ha="center", va='top',
                 arrowprops=dict(shrink=0.05, facecolor='k'))

#凡例の表示。Figure の判例を出力することもできる。
ax1.legend(loc='upper right', frameon=False)

gs.tight_layout(fig) #グラフ描画前に行うべし。慣例では、fig.tight_layout()
plt.show(fig) #慣例では plt.show() と fig を指定しない。

ValueError: operands could not be broadcast together with shapes (1500,) (1500,3)
このエラーが原因だと思うのですが、どうしたらいいかわからないです。詳しい方お力をお借りしたいです。

Comments

[melian]

shape が (1500,3) なので、画像を BGR モードで読み込んでいるのだと思います。なので、グレースケールで読み込むとよいのではないでしょうか。
img = cv2.imread('img.png', 0)

それと、quadratic 関数のパラメータ x にリストが渡される場合があるので、そこでエラーが発生します。

[mappys]

melian様
回答ありがとうございます！
img = cv2.imread('img.png', 0)
で読み込んでみたところ、エラー内容が変わりました。質問させていただいたエラーはクリアできたのかなと思います。ありがとうございます！

次に出たエラーが以下になります。

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/scipy/optimize/minpack.py:833: OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',
Traceback (most recent call last):
File "/Users/masashi/pyworks/gaus.py", line 81, in <module>
ax1.plot(channels, fitting(channels, *paramater_optimal),
File "/Users/masashi/pyworks/gaus.py", line 36, in fitting
return quadratic(x, a0, b0, c0) + \
File "/Users/pyworks/gaus.py", line 23, in quadratic
return a + b * x + c * x**2
TypeError: can't multiply sequence by non-int of type 'numpy.float64'

このようなまた別のエラーが出てしまいました。。。

>>それと、quadratic 関数のパラメータ x にリストが渡される場合があるので、そこでエラーが発生します。
ここに関していまいち理解できていないのですが、ご指摘いただいたこの点がエラーになってしまっている感じでしょうか？

[melian]

はい、そうです。quadratic 関数の中身をみますと、x は numpy の array でないとエラーになってしまいます。なので、

def quadratic(x, a, b, c):
x = np.array(x)
return a + b * x + c * x**2

としてみるとよいかと思います。これでフィッティングの結果はでますが。。。

[mappys]

エラー文を見なくてもエラーが出るところがわかるのがすごすぎて感動しております。
ご指摘いただいた部分を書き加えたところ無事結果が出力されました。
ただまともにフィッティングできていない様子です笑
初期値などのパラメータが適切ではないためでしょうか。。。
結果は出力されますが以下の文が出てきます

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/scipy/optimize/minpack.py:833: OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',

おそらくパラメータが適切ではないよ、って意味ですよね。

コードの初期値の部分↓
initial_quadratic = 0.1, 0.1, 0.1 #a0, b0, c0
initial_gauss1 = 2500, 464, 170. #a1, sigma1, mean1
initial_gauss2 = 2400, 1000, 100. #a2, sigma2, mean2
initial_gauss3 = 2000, 1300, 100. #a3, sigma3, mean3

ここを見直す必要がありそうですね。
ひとまず出力まではいけました！本当にありがとうございます。

[mappys]

出力された結果を見ると、quadratic 関数は特に必要なく、ガウス関数のみでフィッテイングした方が良さそうな気がしてきましたので、quadratic 関数を消す方向でもやってみようと思います。

[mappys]

試行錯誤した結果、無事フィッティングできました！
本当にありがとうございました。

Answer 1

melian様のアドバイスで解決できました。
・画像の読み込みをグレースケールで行うこと
・def quadratic(x, a, b, c):
x = np.array(x)
return a + b * x + c * x**2

と書き足すことで解決しました。