分析ノート

もう結構古い記事なのですが、以前SciPyで定積分をやる方法を記事にしたことがあります。
参考: scipyで定積分

最近、2変数関数の積分をやる機会があったのでこの機会に重積分をSciPyで行う方法を紹介します。SciPyのintegrateモジュールには、重積分用の関数が複数あります。
– dblquad (2変数関数の定積分)
– tplquad (3変数関数の定積分)
– nquad (一般のn変数の定積分)

dblquadの使い方

順番に説明していきます。まずは2重積分のdblquadです。関数の定義は次のようになっています。
scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

必須なのは、積分対象のfunc, 外側の積分区間のa, b、そして内側の積分区間を示す、gfun, hfunです。

gfunとhfunは名前からわかる通り、定数ではなく関数です。これにより内側の積分の積分区間を変数にすることができます。つまり以下のような積分区間の積分ができます。
$$\int_0^1\int_0^y xy \,dxdy$$

例に挙げたのでこれを実装してみましょう。ちなみに解は$1/8=0.125$です。funcの定義は、内側の関数を第1引数にする必要があるので注意してください。

from scipy.integrate import dblquad


def f1(x, y):
    return x*y


def x0(y):
    return 0


def x1(y):
    return y


print(dblquad(f1, 0, 1, x0, x1))
# (0.125, 5.515032205777789e-15)

想定通りですね。積分結果と推定誤差が返ってくるのは1変数の積分と同様です。

内側の積分区間も定数から定数までだよ、要するに長方形領域で積分したいよ、って場合はgfun, hfun に定数を返す関数を返してください。

tplquadの使い方

続いて、3変数向けのtplquadです。これもdblquadとかなり似てる感じで使えます。積分変数が一個増えているので上限加減の指定がもう一個ある感じです。
scipy.integrate.tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

たとえば次の積分をやってみましょう。

$$\int_0^1\int_0^z\int_0^{y+z} xyz \,dxdydz.$$

ちなみに答えは$17/144=0.11805555…$となるはずです。

引数の順番に注意が必要なので慎重にコーディングしてください。

from scipy.integrate import tplquad


def f2(x, y, z):
    return x*y*z


def x0(y, z):
    return 0


def x1(y, z):
    return y + z


def y0(z):
    return 0


def y1(z):
    return z


print(tplquad(f2, 0, 1, y0, y1,  x0, x1))
# (0.11805555555555557, 2.1916761217856673e-14)

バッチリですね。

nquadの使い方

最後に一般のn変数を積分できるnquadの使い方を紹介します。

引数の形式が先ほどの二つと少し違います。
scipy.integrate.nquad(func, ranges, args=None, opts=None, full_output=False)

funcにn変数関数を渡して、rangesに積分区間を渡すことになります。rangesは配列で、1変数目から順番に区間の下限上限の2値の配列を格納しておけば良いです。また、ここにも一応関数を使うことはできます。

これはシンプルな例で、定数関数1を超立方体区間で積分してみました。

from scipy.integrate import nquad


def f3(w, x, y, z):
    return 1

print(nquad(f3, [[-2, 2], [-2, 2], [-2, 2], [-2, 2]]))
# (256.0, 2.8421709430404007e-12)

$4^4=256$になりましたね。

ここで急にシンプルな例を出したのには事情がありまして、変数の数が多くなるとやはり積分は困難なようで、ちょっと複雑な例になると規定の反復回数をこなしても必要な精度に届かずWarningが出たりするケースが多々あります。

どうしても計算したい場合は limit パラメーター等をいじっての対応になりますのでドキュメントを参照しながら調整してみてください。
(僕も実運用で必要になったら改めて調査して紹介しようと思います。)

Pythonの多くのライブラリの様々な関数が非常に汎用的に使えるように作られているので多くの引数を受け取れるようになっています。しかし、そのほとんどの引数を固定して1変数関数として使いたいなぁと思うようなことがあります。PandasのDataFrameのapplyなど関数を引数として受け取る関数に渡す場合等ですね。
また、大量にある引数のほとんどを固定して一部だけ変えながら何度も実行する、といった場面も考えられます。

lambda式などを作ってラップした新しい関数を実装してもいいのですが、 functoolsという標準ライブラリにその専用のpartial というメソッドが用意されています。
参考: functools.partial(func, /, *args, **keywords)

このpartialを使うと、引数の一部を固定した引数の少ない新しい関数を作ってくれます。

一個目の引数に元になる関数を渡し、2個目以降の引数に渡したものが、元の関数の固定引数として使われます。keyword引数で渡せばそのkeyword引数が固定されます。

一引数の固定の方は先頭から順番に固定されるので注意してください。つまり2番目以降の引数を固定したい場合はそれらはキーワード引数として指定する必要があります。

サンプル

引数を順番に表示するだけの単純な関数を作ってやってみましょう。

from functools import partial


# 3つの引数を表示するだけの関数
def sample_func(a, b, c):
    print("a=", a)
    print("b=", b)
    print("c=", c)


# テスト実行
sample_func(1, 2, 3)
"""
a= 1
b= 2
c= 3
"""

# a = 10, b = 20 を固定した新しい関数が作られる。
partial_f = partial(sample_func, 10, 20)


# 3個目の引数 c = 50だけ渡して実行できる。
partial_f(50)
"""
a= 10
b= 20
c= 50
"""

# キーワード引数で固定することもできる。
partial_f2 = partial(sample_func, a=100, c=200)

# b の値だけ渡して実行できる
partial_f2(b=-5)
"""
a= 100
b= -5
c= 200
"""

キーワード引数を固定した関数を、位置引数で使う場合は注意が必要です。
たとえば、次のようにaを固定して生成した関数に、残り2個の引数を位置引数で渡すと、aを2回渡した扱いになってエラーが起きます。

# aを固定
partial_f3 = partial(sample_func, a=1)

# bとcのつもりで残り2個の引数を渡すとエラー
try:
    partial_f3(2, 3)
except Exception as e:
    print(e)
# sample_func() got multiple values for argument 'a'

# bとcもキーワード引数で渡す。
partial_f3(b=2, c=3)
"""
a= 1
b= 2
c= 3
"""

まとめ

ほぼ小ネタのような内容でしたが、自作関数をベースに一部の振る舞いを固定した簡易的な関数を作るとか、apply等の1変数関数を受け取るメソッドに渡したいとかそういう場面で役に立つことがあるテクニックとしてpartialを紹介しました。

scipyのstats配下の各種メソッドであれば、それぞれがパラメーターを固定するfrozenメソッドを持ってるとか、引数が多いなら引数を辞書にまとめて**(アスタリスク2個)で展開すればいいとか、ラップした関数を自分で実装したらいいとか、代用手段も多いのですが、partialを使うとその辺の記述がシンプルになるので機会があれば使ってみてください。

前回の記事の二分法に続いて、もう一つ求根アルゴリズムを紹介します。
参考: 二分法を用いて関数の根を求める

今回紹介するのはニュートン法です。これは微分可能な関数$f(x)$の根を求めることができるアルゴリズムです。
参考: ニュートン法 – Wikipedia

詳しい説明は上記のWikipediaにあるので、ざっくりと概要を説明します。

この方法の背景にあるのは、滑らかな関数をある点の近くだけ着目してみるとほぼ直線になり、接線で近似できるということをベースのアイデアにしています。

つまり、微分可能な関数$f$があって$f(x)=0$だとします。その根$x$の近くに点$x_0$を取ると、$x_0$の近くでは、$f$と$f$の接線ってかなり近いよね、それなら$f$の根と$f$の$x_0$における接線の根って近いよね、っていうのが基本的なアイデアです。

関数$f$の$f(x_0)$における接戦は次の式で書けます。

$$y=f'(x_0)(x-x_0)+f(x_0).$$

$f(x)=0$は解けない場合でも、この接線の根は容易に算出することができ、

$$x_0-\frac{f(x_0)}{f'(x_0)}$$

と求まります。

この値は元の$x_0$よりも真の根に近いことが期待され、これをもう一回$x_0$とおいて同じ操作を繰り返せば真の根にたどり着く、というのがニュートン法です。

Wikipediaから画像拝借しますが、図で見るとイメージしやすいですね。

注意しないといけないのは、初期値$x_0$は真の根$x$の十分近くに取らないといけない点です。十分近くを見れば関数をその接線で近侍できるよね、というのがアイデアの前提なので、根が近くになかったらその前提が崩れてしまいこのアルゴリズムは真の根に収束しなくなってしまいます。

ニュートン法のメリットとデメリット

先に紹介した二分法と比べて、ニュートン法のメリットデメリットを説明していきます。

１番のメリットは収束の速さです。二分法に比べてより少ない計算回数で効率的に会を探索することができます。

また、初期値として与える点が1点だけで良いというのもメリットです。二分法の場合は初期値は区間で設定する必要がありましたからね。

その一方で複数のデメリットもあります。実装していて一番不便に感じるのはその関数だけでなく微分も必要ということでしょうか。もちろん微分不可能な関数ではニュートン法は使えません。

また、初期値が真の解に十分近くない場合や、微分した値が$0$に近い場合、うまく収束せずにアルゴリズムが失敗してしまう、という点も大きなデメリットです。

SciPyによる実装

SciPyではscipy.optimizeというモジュールで実装されています。newtonという専用メソッドを使うか、root_scalarという汎用的なメソッドで(method=’newton’)を指定して使うことになります。二分法と同じですね。

参考:
scipy.optimize.newton — SciPy v1.12.0 Manual
root_scalar(method=’newton’) — SciPy v1.12.0 Manual

二分法の時と同じように、$\sin$関数の根$\pi$を探索させてみましょう。微分は$\cos$なのでこれを使います。

from scipy import optimize
import numpy as np


root1 = optimize.newton(np.sin, x0=3, fprime=np.cos)
print(root1)
# 3.141592653589793

root_result = optimize.root_scalar(np.sin, method="newton", x0=3, fprime=np.cos)
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 6
     iterations: 3
           root: 3.141592653589793
"""

print(root_result.root)
# 3.141592653589793

簡単ですね。

注目するのは、iterationsの部分です。たった3回のイテレーションで収束していて、関数が実行されたのは、fとfの微分合わせて6回だけです。
二分法の時は39回もイテレーションが必要だったのと大違いです。そして実はこの例では解の精度もニュートン法の方が高くなっています。

ニュートン法が失敗する例

初期値が真の解の近くにないと失敗するという話がありましたのでそちらも見ておきます。

例えば、タンジェントの逆関数、$\arctan$で試してみましょう。(sin, cosは根が無限にあって、根から遠い実数を用意できないので関数を変えます。)

$\arctan(x)$の微分は$\frac{1}{1+x^2}$です。

やってみました。

def f(x):
    return np.arctan(x)


def fprime(x):
    return 1/(1+x**2)


# 初期値が1なら収束する
root_result_1 = optimize.root_scalar(np.sin, method="newton", x0=1, fprime=fprime)
print(root_result_1)
"""
      converged: True
           flag: 'converged'
 function_calls: 12
     iterations: 6
           root: 0.0
"""

# 初期値が2だと失敗し、結果のflagが'convergence error'になる。
root_result_2 = optimize.root_scalar(np.sin, method="newton", x0=2, fprime=fprime)
print(root_result_2)
"""
      converged: False
           flag: 'convergence error'
 function_calls: 100
     iterations: 50
           root: 1.854706857103781
"""


# optimize.newton の方だと例外が上がる。
try:
    optimize.newton(f, fprime=fprime, x0=2)
except Exception as e:
    print(e)
# Derivative was zero. Failed to converge after 10 iterations, value is -6.999943395317963e+168.

失敗した時の振る舞いがそれぞれ違うので、どちらのコードを使うかで注意深く扱う必要がありますね。optimize.root_scalarはコード自体は正常に終了しますがフラグが立ち、optimize.newtonの方は例外があがります。

ちなみに、エラーの中で出てくるvalue の値、 -6.999943395317963e+168 は次のように自分でニュートン法を実装しても同じ値が出て来ます。

x0 = 2  # 初期値
for i in range(12):
    x0 = x0 - f(x0)/fprime(x0)
    print(i+1, "回目: x0=", x0)

"""
1 回目: x0= -3.535743588970453
2 回目: x0= 13.950959086927496
3 回目: x0= -279.34406653361754
4 回目: x0= 122016.9989179547
5 回目: x0= -23386004197.933937
6 回目: x0= 8.590766671950415e+20
7 回目: x0= -1.1592676698907411e+42
8 回目: x0= 2.110995587611039e+84
9 回目: x0= -6.999943395317963e+168
10 回目: x0= inf
11 回目: x0= nan
12 回目: x0= nan
"""

絶対値が大きくなり続けていて全く収束に向かっていないのがわかりますね。

まとめ

元の関数だけではなく導関数も必要だったり、初期値の設定段階である程度解の目星をつけておかないといけないなどのデメリットはありますが、速度や精度の面で優秀でしかもロジックもわかりやすい手法なので、何か機会があればニュートン法の活用を検討してみてください。

1変数連続関数の根(値が0になる点)を求める、二分法というアルゴリズムとそれをScipyで実装する方法を紹介します。

アルゴリズムの内容

二分法というのは中間値の定理をベースとした求根アルゴリズムです。アイデアは非常に単純で、連続関数$f$に対して、$f(x_1)$と$f(x_2)$の符号が異なるように、$x_1, x_2$を選びます。この時点で、中間値の定理より区間$(x_1, x_2)$に根があることがわかりますのでさらに細かくみていきます。次は、$x_1, x_2$の中点$x_M = \frac{x_1+x_2}{2}$を取り、$f(x_M)$の符号を調べます。$f(x_M)$と$f(x_1)$の符号が同じであれば、$x_1$を$x_M$で置き換え、逆に$f(x_M)$と$f(x_2)$の符号が同じであれば、$x_2$を$x_M$で置き換えると、区間の幅が半分になった$(x_1, x_2)$が得られますが根はこの中にあることがわかります。これを繰り返すと、根が存在する範囲を狭めていくことができ、$f(x_M)$の絶対値が0になるか、もしくは十分0に近づいたらその値を数値的に求めた根とします。

以上が、一般的な二分法のアルゴリズムの説明です。ただし、後に紹介するSciPyではどうやら区間が十分狭くなったかイレーション回数が上限に達したか等の基準でループを打ち切り、$f(x)$の値を確認していないようです。

二分法のメリットとデメリット

方法が単純でわかりやすい、というのが個人的に感じている１番のメリットです。

また、連続関数であれば使えるため、ニュートン法などのアルゴリズムのように元の関数の微分を必要とせず、微分が難しい関数や微分不可能な関数でも使えます。

また、根が存在しうる区間を狭めながら探索するため、最初の区間の幅と繰り返し回数により、結果の精度を保証できることも大きな利点です。

逆にデメリットとしては、ニュートン法等と比較して収束が遅いこととか、初期値として関数が異符号になる2点を探して与える必要があること、もしその区間に複数の根が存在した場合にどれに収束するか不確定なことなどが挙げられます。

ただし僕の経験上では、ある程度根の目処がついていたり、単調な関数に対して使う場面が多くこれらのデメリットを深刻に感じることは少ないです。

SciPyによる実装

SciPyではscipy.optimizeというモジュールで実装されています。bisectという専用メソッドを使うか、root_scalarという汎用的なメソッドで(method=’bisect’)を指定して使うことになります。

参考:
– scipy.optimize.bisect — SciPy v1.11.4 Manual
– root_scalar(method=’bisect’) — SciPy v1.11.4 Manual

試しに、$\sin$関数の 3と4の間にある根を探させて見せましょう。既知の通りそれは円周率$\pi$になるはずです。

from scipy import optimize
import numpy as np


root1 = optimize.bisect(np.sin, 3, 4)
print(root1)
# 3.1415926535901235

root_result = optimize.root_scalar(np.sin, bracket=[3, 4], method='bisect')
# 結果は複数の情報を含むRootResults形式で戻る。
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 3.1415926535901235
"""
print(type(root_result))
# class 'scipy.optimize._zeros_py.RootResults'

# 根の値へのアクセス方法
print(root_result.root)
# 3.1415926535901235

いかにも円周率ぽい結果が得られましたね。root_scalarの方では収束したことを示すフラグや、イテレーション回数なども得られています。

失敗事例1. 区間の両端での関数の値が同符号の場合

二分法は初期設定を誤ってると失敗するのでその場合のSciPyの挙動も見ておきましょう。失敗パターンの一つは、最初に指定した区間の両端で符号が一致していた場合です。もちろん関数の形によってはその区間内に根がある可能性もあるのですが、存在は保証されなくなります。

また$\sin$関数で、その間に根が存在しない区間$(1, 2)$と、実は両端で同符号だけど根が存在する区間$(1, 7)$でやってみましょう。bisectとroot_scalarで全く同じエラーメッセージ出るのでbisectの方だけ載せます。

try:
    optimize.bisect(np.sin, 1, 2)
except Exception as e:
    print(e)
# f(a) and f(b) must have different signs

try:
    optimize.bisect(np.sin, 1, 7)
except Exception as e:
    print(e)
# f(a) and f(b) must have different signs

はい、$f(a)$と$f(b)$は違う符号にせよとのことでした。根が存在しない1個目の例はさておき、2個目の例は根は区間内に2個存在するのですが探さずにエラーになりました。

失敗事例2. 連続関数ではなかった場合

もう一つ失敗するのは、関数が連続関数ではないケースです。

例えば$\tan$の、$\frac{\pi}{2}$近辺の挙動で見てみましょう。数学的に厳密な話をすると、$\frac{\pi}{2}$では$\tan$は定義されないので、$\tan$は数学的には連続関数(定義域内のすべての点で連続)なのですが、数値計算的には不連続と考えた方が都合が良いです。

話が脇に逸れたので実例の話に移ります。実は、bisectメソッドは結果を返して来ちゃうんですよね。そしてそれが全然根ではないということも見ておきましょう。

root = optimize.bisect(np.tan, 1, 2)
# pi/2に近い結果が得られている
print(root)
# 1.5707963267941523

# 元の関数に代入した結果は全く0に近くない
np.tan(root)
# 1343445450736.3804

root_scalarの方だったら、結果のフラグ等もあるのでアラート等あげてくれるのかと期待したのですがそういう機能はなさそうです。

root_result = optimize.root_scalar(np.tan, bracket=[1, 2], method='bisect')
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 1.5707963267941523
"""
np.tan(root_result.root)
# 1343445450736.3804

要するに、SciPyに渡す関数が本当に連続関数であるかどうかは利用者が責任を持たないといけないということです。また、結果が本当に根なのかどうかは代入して確認した方が良いでしょう。

まとめ

以上が二分法の説明とSciPyで利用する方法、その注意点でした。