質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.37%

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

解決済

回答 3

投稿

  • 評価
  • クリップ 0
  • VIEW 1,134

throwsNullPo

score 12

解決したいこと

変更禁止の関数(下記の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))
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 3

+3

np.frompyfuncnp.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))

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2019/06/10 22:25

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

    キャンセル

checkベストアンサー

+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]]))

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2019/06/07 10:06

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

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

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

    キャンセル

  • 2019/06/10 22:25

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

    キャンセル

0

単に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()

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2019/06/07 04:15

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

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.37%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る