scipyのcurve_fitを変更禁止の関数に対して適用したい

Question

解決したいこと

変更禁止の関数（下記のtarget_func2）を直接的もしくは間接的にcurve_fitに渡して推定結果を得たいのですが、test_fit2を実行すると以下のエラーが出ます。
ValueError: callable <ufunc '? (vectorized)'> is not supported by signature
target_func2を一切変更せずにtest_fit1のような結果を得るにはどうすればいいでしょうか。

サンプルデータはa=3.0にして、そこからちょっとずれるように作っています。test_fit1の結果は
(array([ 3.00033445]), array([[ 3.35566718e-05]]))
で想定通りでした。

from scipy import optimize
import numpy as np


def target_func1(x, a):
    import scipy
    return (scipy.exp(x) - 1.0)/x + a

def test_fit1():
    xdata = np.arange(*[0.01, 3.0, 0.01]) #とりあえず動かすため0を含めない
    ydata = np.array([target_func1(x, 3.0) + (-1.0)**i / 10.0 for i, x in enumerate(xdata)])
    print(optimize.curve_fit(target_func1, xdata, ydata, maxfev=1000))


def target_func2(x, a):
    if abs(x) < 1.0e-6:
        return 1 + x/2 + x**2/6 + x**3/24 + x**4/120 + a
    
    import math
    return (math.exp(x) - 1.0)/x + a

def test_fit2():
    xdata = np.arange(*[0.0, 3.0, 0.01]) # 0を含む
    ydata = np.array([target_func2(x, 3.0) + (-1.0)**i / 10.0 for i, x in enumerate(xdata)])
    univ_func = np.frompyfunc(target_func2, 2, 1)
    print(optimize.curve_fit(univ_func, xdata, ydata, maxfev=1000))

Answer 1

np.frompyfuncをnp.vectorizeで置き換えると動くようです。
置き換えただけでは、ValueError: Unable to determine number of fit parameters.がでるので、curve_fitのオプションにp0=[1](パラメータの初期値)を追加しています。

def test_fit2():
    xdata = np.arange(*[0.0, 3.0, 0.01]) # 0を含む
    ydata = np.array([target_func2(x, 3.0) + (-1.0)**i / 10.0 for i, x in enumerate(xdata)])
    univ_func = np.vectorize(target_func2, excluded=[1])
    print(optimize.curve_fit(univ_func, xdata, ydata, p0=[1], maxfev=1000))

Comments

[throwsNullPo]

ありがとうございます。無事意図したとおりに動かすことができました。

Answer 2

curve_fitはターゲットの関数オブジェクトを調べそこから引数情報を引き出してfit parametersをどのように制御するかを決める仕組みになっているようです。ところがnumpyのuniversal functionにしてしまうと元の関数の引数情報が失われ、ご質問のエラーとなるっているようでした。ラッパー関数によくある「*args, **kwargs」という汎用的な仮引数の記述がオリジナルの関数の引数情報を消してしまっているということだと思います。

つまりnumpy.frompyfuncで生成したユニバーサル関数はcurve_fitにとって都合の悪い実装でありそれがhayataka2049さん回答にある「numpyのユニバーサル関数にしているのが駄目」の意味なのだと思いました。そこで型がcurve_fitにわかるようにラッパー関数を書いてみました。

前述の問題があるため普通のラッパー関数でするように*args, **kwargsといった仮引数の書き方はできないためオリジナル関数の引数パターンをそのまま踏襲した特別仕立てのラッパーにしたわけです。

...

def test_fit2():
    xdata = np.arange(*[0.0, 3.0, 0.01]) # 0を含む
    ydata = np.array([target_func2(x, 3.0) + (-1.0)**i / 10.0 for i, x in enumerate(xdata)])
    # univ_func = np.frompyfunc(target_func2, 2, 1)
    univ_func = ufunc_with_arg2(target_func2)
    print(optimize.curve_fit(univ_func, xdata, ydata, maxfev=1000))


def ufunc_with_arg2(f):
    ufunc_impl = np.frompyfunc(f, 2, 1)

    def ufunc(arg1, arg2):
        return ufunc_impl(arg1, arg2)

    return ufunc

これでうまくいくと思いきや・・・

minpack.error: Result from function call is not a proper array of floats.

と言われしまいました。あれれ・・・と思って上記のufunc_implの結果の型を調べてみるとndarrayにはなっているのですがdtypeがobjectになってました。え？っと思ってリファレンスを見ると

https://docs.scipy.org/doc/numpy/reference/generated/numpy.frompyfunc.html

Notes:
The returned ufunc always returns PyObject arrays.

なるほどこいつは数値に限らず任意のデータ型用のユニバーサル関数を作るためのもののようです。そこでnumpy.frompyfuncじゃなくてnumpy.vectorizeを使ってみると今度はエラーなく結果がでました。こちらは数値演算が前提のブロードキャスティング対応のラッパー関数を作ってくれるということだと思います。

...

def ufunc_with_arg2(f):
    # ufunc_impl = np.frompyfunc(f, 2, 1)
    ufunc_impl = np.vectorize(f)

    def ufunc(arg1, arg2):
        return ufunc_impl(arg1, arg2)

    return ufunc

...

==> (array([3.]), array([[3.34448161e-05]]))

これでいかがでしょうか？

Comments

[KSwordOfHaste]

ちなみに自分はscipyの機能をほとんど何も知らないので、もっと平易な方法があるのではないかとも想像しました。つまり

curve_fit_ex(func, num_parameters, ...)のようにしてfuncの第一引数以外のパラメータの数を指定できる・・・

といったようなものです。 調べてないのでそういうものがあるかどうかは分かりませんが・・・

[throwsNullPo]

詳しい説明ありがとうございます。落ちている個所までデバッグしたものの何をやっているのか最初はわからなかったのですがおかげで腑に落ちました。

Answer 3

単にnumpyのユニバーサル関数にしているのが駄目なので、これだけ外せればうまくいきます。

pythonのループに書き換えましょう。下のコードは内包表記の例です。

import math
import numpy as np
from scipy import optimize

def target_func2(x, a):
    if abs(x) < 1.0e-6:
        return 1 + x/2 + x**2/6 + x**3/24 + x**4/120 + a

    return (math.exp(x) - 1.0)/x + a

def test_fit2():
    xdata = np.arange(*[0.0, 3.0, 0.01])
    ydata = np.array([target_func2(x, 3.0) + (-1.0)**i / 10.0
                      for i, x in enumerate(xdata)])

    result = [optimize.curve_fit(target_func2, x, y, maxfev=1000) 
              for x, y in zip(xdata, ydata)]

test_fit2()

Comments

[throwsNullPo]

回答ありがとうございます。
上記コードを張り付けて実行してみましたが、思っていたのと違った結果が返ってきました。
上の場合、サンプルから1ペアずつ取り出して、シングルデータでフィッティングを行っているように見えます。最初の要素は
optimize.curve_fit(target_func2, 0.0, 4.0+0.1, maxfev=1000)
と等価と思われ、実際この結果は
(array([ 3.1]), array([[ inf]]))
となり、張っていただいたコードのresultの最初の要素と一致しました。