2次元ガウス関数をdblquadで積分したものでfittingを行いたい。ValueError: The truth value of an array with more than ...

Question

以下のデータを二次元ガウス関数の積分を用いてfittingをしたいと考えています。このデータはある観測データであり、5×5の二次元データを一次元に均したものです。主な構造として5×5の中心の値がほかの値より大きくなっており、ここでは2738.33がそれにあたります。

z = [359.44, 310.44, 317.78, 314.78, 323.0, 316.44, 416.56, 330.0, 322.33, 363.89, 290.11, 312.22, 2738.33, 343.89, 323.89, 336.0, 355.67, 359.44, 327.67, 336.33, 319.67, 332.0, 277.0, 386.56, 344.33]

二次元ガウス関数の積分を以下のようにします。
各パラメータの説明をします。heightはガウス関数にかける数で、ガウス関数の体積は1となるので、このheightの値はガウス関数の体積を表しています。sigmaはガウス関数のなだらかさを示すパラメータ。cen_xとcen_yはガウス関数の頂点が中心座標の(0, 0)からどれだけ離れているかを表します。constantは初めから与えられる初期値を表します。これらのパラメータをfittingで定めていきます。
積分はintegrate.dblquadを用いて行っていきます。x_from, x_to, y_from, y_toにはそれぞれ積分範囲のxの始まりと終わり、yの始まりと終わりを表します。integrate.dblquadは（積分値、推定誤差）を出力するのでここで定義した関数では第一引数のみ出力します。

def quad_gaussian2D(param, x_from, x_to, y_from, y_to):
    height = param[0]
    sigma = param[1]
    cen_x = param[2]
    cen_y = param[3]
    constant = param[4]
    quad_result = integrate.dblquad(lambda x, y: height/(2*np.pi*sigma**2)*np.exp(-(((cen_x-x)**2+(cen_y-y)**2)/(2*sigma**2))) + constant, x_from, x_to, y_from, y_to)
    return quad_result[0]

積分範囲の設定を行います。
以下のように長さが25の４つのリストを作成し、それぞれがx_from, x_to, y_from, y_toに対応します。
quad_range[0]
[-2.5 -2.5 -2.5 -2.5 -2.5 -1.5 -1.5 -1.5 -1.5 -1.5 -0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 0.5 0.5 0.5 1.5 1.5 1.5 1.5 1.5]
quad_range[1]
[-1.5 -1.5 -1.5 -1.5 -1.5 -0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 0.5 0.5 0.5 1.5 1.5 1.5 1.5 1.5 2.5 2.5 2.5 2.5 2.5]
quad_range[2]
[-2.5 -1.5 -0.5 0.5 1.5 -2.5 -1.5 -0.5 0.5 1.5 -2.5 -1.5 -0.5 0.5 1.5 -2.5 -1.5 -0.5 0.5 1.5 -2.5 -1.5 -0.5 0.5 1.5]
quad_range[3]
[-1.5 -0.5 0.5 1.5 2.5 -1.5 -0.5 0.5 1.5 2.5 -1.5 -0.5 0.5 1.5 2.5 -1.5 -0.5 0.5 1.5 2.5 -1.5 -0.5 0.5 1.5 2.5]

q_list = []
for i in np.arange(-2.5, 1.6):
    q_list.append([i, i+1])
quad_range = []
for i in q_list:
    for j in q_list:
        quad_range.append(i+j)
quad_range = np.array(quad_range)
quad_range = quad_range.T

x_from = quad_range[0]
x_to = quad_range[1]
y_from = quad_range[2]
y_to = quad_range[3]

fittingをおこなっていきます。
まずパラメータの初期値を決めていきます。
heightはデータの中心の値、constantはデータの最外のデータの平均値を採用しました。
residuals_2Dでデータとガウス関数の積分値の引き算を行う。今回は積分区画がx, yで範囲1ずつであるので、観測データの値を高さとした長方体を考えたときにデータ値がそのまま体積となるので、ガウス関数の積分値との引き算を行っている。fittingではこの値が最も小さくなるパラメータを探す最小二乗法で行っている。

param_2D = [2738.33, 1, 0, 0, 328.22875]

def residuals_2D(param, z, x_from, x_to, y_from, y_to):
    return (z - quad_gaussian2D(param, x_from, x_to, y_from, y_to))

optimised_param_2D = optimize.leastsq(residuals_2D, param_2D, args = (z, x_from, x_to, y_from, y_to))
height_2D = optimised_param_2D[0][0]
sigma_2D = optimised_param_2D[0][1]
cenx_2D = optimised_param_2D[0][2]
ceny_2D = optimised_param_2D[0][3]
cons_2D = optimised_param_2D[0][4]

コード全容

import numpy as np
from scipy import optimize
from scipy import integrate

def quad_gaussian2D(param, x_from, x_to, y_from, y_to):
    height = param[0]
    sigma = param[1]
    cen_x = param[2]
    cen_y = param[3]
    constant = param[4]
    quad_result = integrate.dblquad(lambda x, y: height/(2*np.pi*sigma**2)*np.exp(-(((cen_x-x)**2+(cen_y-y)**2)/(2*sigma**2))) + constant, x_from, x_to, y_from, y_to)
    return quad_result[0]

def residuals_2D(param, z, x_from, x_to, y_from, y_to):
    return (z - quad_gaussian2D(param, x_from, x_to, y_from, y_to))

z = [359.44, 310.44, 317.78, 314.78, 323.0, 316.44, 416.56, 330.0, 322.33, 363.89, 290.11, 312.22, 2738.33, 343.89, 323.89, 336.0, 355.67, 359.44, 327.67, 336.33, 319.67, 332.0, 277.0, 386.56, 344.33]

q_list = []
for i in np.arange(-2.5, 1.6):
    q_list.append([i, i+1])
quad_range = []
for i in q_list:
    for j in q_list:
        quad_range.append(i+j)
quad_range = np.array(quad_range)
quad_range = quad_range.T

x_from = quad_range[0]
x_to = quad_range[1]
y_from = quad_range[2]
y_to = quad_range[3]

param_2D = [2738.33, 1, 0, 0, 328.22875]

optimised_param_2D = optimize.leastsq(residuals_2D, param_2D, args = (z, x_from, x_to, y_from, y_to))
height_2D = optimised_param_2D[0][0]
sigma_2D = optimised_param_2D[0][1]
cenx_2D = optimised_param_2D[0][2]
ceny_2D = optimised_param_2D[0][3]
cons_2D = optimised_param_2D[0][4]

このコードを実行したときに
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
が発生してしまいます。デバッグを行うとquad_gaussian2D関数のquad_result =の行で発生してることがわかりました。
このエラーはbool値を要素とするnumpyを使用したときに発生するエラーであると認識しているのですが、自分ではbool値による処理を行っているつもりはないので、どこが間違っているのかわかりません。どこを直せば改善するかお分かりになる方がいらっしゃったらご指摘していただけると幸いです。

Answer 1

エラーが起こっている原因は、 integrate.dblquad に配列を渡しているからです。
dblquadの始点・終点はスカラーを渡しますが、配列を渡しているため、中で比較したときにそのエラーが発生しています。

residuals_2D() で引数をprintしたらわかると思いますが、z, x_from, x_to, y_from, y_to 等は、ひとつひとつスカラーで呼び出されるのではなく、配列が一回で渡されます。
その配列をそのままquad_gaussian2D()に渡して、integrate.dblquadを呼び出しているのが問題です。

解決策としては、例えば下記のように、forでひとつひとつquad_gaussian2D()を実行すればいいです。
(速度的にはいまいちなやり方ですが)

def residuals_2D(param, z, x_from, x_to, y_from, y_to):
    # print(param, x_from, x_to, y_from, y_to)
    return [
        z0 - quad_gaussian2D(param, x0, x1, y0, y1)
        for z0, x0, x1, y0, y1 in zip(z, x_from, x_to, y_from, y_to)]

Comments

[python_wakaran]

回答に気づかず申し訳ございません。
以上のように改善したところ、解決いたしました。
ありがとうございました。