6.17. ソルバーと乱数

対象とする関数が、サンプルのバッファではなくPythonコードによって定義されている場合は、別系統のツールが適切です。関数の根はどこか、その最小値は、与えられた区間にわたる積分は。scipy.integratescipy.optimize のサブモジュールがこの作業をカバーします。各アルゴリズムはユーザーが提供したPython関数にコールバックするため、反復ごとのコストはバッファのリダクションよりも高くなります。利点は、ソルバーを自分で書かなくてよいことです。

scipy.special サブモジュールは、確率分布の累積分布関数(CDF、サンプルが特定の値以下である確率)や確率密度関数(PDF、特定の値における相対的な尤度)を計算する際に出てくる統計的特殊関数(誤差関数、ガンマ関数)をカバーします。numpy.random は、ディザリング、シミュレーション、合成テストデータ用の疑似乱数生成器をカバーします。

6.17.1. 呼び出し可能オブジェクトの数値積分

被積分関数がサンプルのバッファではなくPython関数である場合、scipy.integrate は4つの求積アルゴリズムを公開します:

  • quad() -- 適応型のガウス・クロンロッド法。滑らかな被積分関数に対する適切なデフォルトです。(value, error) を返します。

  • romberg() -- 古典的なロンバーグ法/ニュートン・コーツ法。単一のfloatを返します。上流では非推奨ですが、互換性のために含まれています。

  • simpson() -- 適応型のシンプソン則。単一のfloatを返します。

  • tanhsinh() -- 二重指数関数型求積法。被積分関数が端点に特異点を持つ場合や、無限の積分限界を持つ場合に使用します。(value, error) を返します。

二重指数関数則(tanhsinh)で評価したガウス積分:

from math import exp
from math import pi
from math import sqrt
from ulab import numpy as np
from ulab import scipy as sp

f = lambda x: exp(-x * x)
value, err = sp.integrate.tanhsinh(f, -np.inf, np.inf)
print("approx:", value, "   exact:", sqrt(pi))

出力:

approx: 1.7724538...   exact: 1.7724538...

6.17.2. 求根と最小化

scipy.optimize は、3つの古典的な単一変数ソルバーをカバーします。各反復はユーザーが提供したPython関数にコールバックするため、純粋なPythonソルバーに対する高速化は控えめ(おおよそ2倍)です。利点は、ソルバーを自分で書かなくてよいことです。

  • bisect() -- 区間を半分にすることで [a, b] 上の f の根を見つけます。f(a)f(b) は符号が逆でなければなりません:

    def f(x):
        return x * x - 1
    
    sp.optimize.bisect(f, 0, 4)        # ~1.0
    
  • newton() -- 割線法/ニュートン・ラフソン反復を使って根を見つけます:

    def f(x):
        return x * x * x - 2.0
    
    sp.optimize.newton(f, 3., tol=0.001, rtol=0.01)
    # ~1.260
    
  • fmin() -- 滑降シンプレックス法(ネルダー・ミード法)を使って局所最小値を見つけます:

    def f(x):
        return (x - 1) ** 2 - 1
    
    sp.optimize.fmin(f, 3.0)           # ~1.0
    

単一変数の範囲で、カメラ側のほとんどの最適化には十分です。センサーのキャリブレーション定数、コントラスト測定を最大化するゲイン、ヒストグラムの双峰性が最も鋭くなるしきい値などです。多変数の問題に対しては、汎用の非線形オプティマイザに頼るよりも、問題を小さな線形代数の求解として再定式化するのが通常は適切な答えです。

6.17.3. 特殊関数

scipy.special は、ユニバーサル関数のように振る舞ういくつかの統計・確率関数を公開します。これらはスカラー、反復可能オブジェクト、または ndarray を受け取り、float ndarray を返します:

x = np.linspace(0, 4, num=8)

sp.special.erf(x)         # error function
sp.special.erfc(x)        # complementary error function
sp.special.gamma(x + 1)   # gamma function
sp.special.gammaln(x + 1) # log-gamma function

誤差関数とその補関数は、ガウス分布のCDFに現れます。測定したzスコアと確率の間を変換したり、正規分布の裾の積分を計算したりするのに最適な用途です。ガンマ関数と対数ガンマ関数は、ベータ/カイ二乗/スチューデントtの計算に登場します。gammaln は、gamma 自体がオーバーフローしてしまう大きな引数に対して数値的に安定な形式です。

6.17.4. 乱数

numpy.random は、一般的な分布からサンプルを引く Generator クラスを提供します。この生成器はステートフルです。各呼び出しが内部状態を進めるため、連続する呼び出しは独立したサンプルを返します:

from ulab import numpy as np

rng = np.random.Generator(seed=42)

rng.random(size=5)             # 5 uniform [0.0, 1.0) samples
rng.uniform(low=-1.0, high=1.0, size=10)
rng.normal(loc=0.0, scale=1.0, size=(2, 4))

出力のdtypeは常に float です。size= は整数(1次元出力)またはタプル(n次元出力)を受け付けます。省略した場合は、単一のPython floatが返されます。

この生成器は、シミュレーション、ディザリング、合成テストデータ、および暗号学的強度が必要とされないその他のあらゆる用途に適しています。鍵やトークンには適していません。それらには os を通じてシステムの乱数源を使用してください。

6.17.5. ビルド時の利用可否

各サブモジュールが実際に存在するかどうかは、カメラがどのようにビルドされたかに依存します。scipy.optimizescipy.special はすべてのカメラで有効になっているわけではありません。カメラに含まれていない関数を呼び出すと AttributeError が発生します。dir(sp)dir(sp.optimize)dir(np.random) などは、対象のカメラで何が利用可能かを報告します。