6.17. Solwery i liczby losowe

Gdy badana funkcja jest zdefiniowana kodem Pythona, a nie buforem próbek, właściwym wyborem jest inna rodzina narzędzi: gdzie leży miejsce zerowe funkcji, jej minimum, jej całka w danym przedziale? Submoduły scipy.integrate i scipy.optimize pokrywają tę pracę. Każdy algorytm wywołuje zwrotnie funkcję Pythona dostarczoną przez użytkownika, więc koszt na iterację jest wyższy niż przy redukcji bufora; wygoda polega na tym, że nie trzeba samemu pisać solwera.

Submoduł scipy.special pokrywa statystyczne funkcje specjalne (funkcja błędu, gamma), które pojawiają się przy obliczaniu dystrybuanty rozkładu prawdopodobieństwa (CDF, prawdopodobieństwo, że próbka jest co najwyżej równa zadanej wartości) lub funkcji gęstości prawdopodobieństwa (PDF, względna wiarygodność przy zadanej wartości). numpy.random pokrywa generator pseudolosowy do ditheringu, symulacji i syntetycznych danych testowych.

6.17.1. Całkowanie numeryczne obiektu wywoływalnego

Gdy funkcja podcałkowa jest funkcją Pythona, a nie buforem próbek, scipy.integrate udostępnia cztery algorytmy kwadratury:

  • quad() – adaptacyjna metoda Gaussa-Kronroda. Właściwy domyślny wybór dla gładkich funkcji podcałkowych. Zwraca (value, error).

  • romberg() – klasyczna metoda Romberga / Newtona-Cotesa. Zwraca pojedynczą liczbę float. Wycofana w źródle nadrzędnym; dołączona dla zgodności.

  • simpson() – adaptacyjna reguła Simpsona. Zwraca pojedynczą liczbę float.

  • tanhsinh() – kwadratura podwójnie wykładnicza. Stosuj, gdy funkcja podcałkowa ma osobliwości na końcach przedziału lub granicę nieskończoną. Zwraca (value, error).

Całka Gaussa obliczona regułą podwójnie wykładniczą (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))

Wynik:

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

6.17.2. Wyznaczanie miejsc zerowych i minimalizacja

scipy.optimize pokrywa trzy klasyczne solwery jednej zmiennej. Każda iteracja wywołuje zwrotnie funkcję Pythona dostarczoną przez użytkownika, więc przyspieszenie względem solwera napisanego w czystym Pythonie jest umiarkowane (mniej więcej 2x); wygoda polega na tym, że nie trzeba samemu pisać solwera.

  • bisect() – znajduje miejsce zerowe f na [a, b] poprzez połowienie przedziału. f(a) i f(b) muszą mieć przeciwne znaki:

    def f(x):
    return x * x - 1

sp.optimize.bisect(f, 0, 4)        # ~1.0

  • newton() – znajduje miejsce zerowe metodą siecznych / Newtona-Raphsona:

    def f(x):
    return x * x * x - 2.0

sp.optimize.newton(f, 3., tol=0.001, rtol=0.01)
# ~1.260

  • fmin() – znajduje minimum lokalne metodą sympleksu schodzącego (Neldera-Meada):

    def f(x):
    return (x - 1) ** 2 - 1

sp.optimize.fmin(f, 3.0)           # ~1.0

Zakres jednej zmiennej wystarcza dla większości optymalizacji po stronie kamery – stała kalibracyjna sensora, wzmocnienie maksymalizujące pomiar kontrastu, próg, przy którym dwumodalność histogramu jest najostrzejsza. W przypadku problemów wielowymiarowych właściwą odpowiedzią jest zazwyczaj przeformułowanie problemu jako małego zadania algebry liniowej, zamiast sięgać po ogólny optymalizator nieliniowy.

6.17.3. Funkcje specjalne

scipy.special udostępnia kilka funkcji statystycznych i prawdopodobieństwa, które zachowują się jak funkcje uniwersalne – przyjmują skalar, obiekt iterowalny lub ndarray i zwracają tablicę 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

Funkcja błędu i jej dopełnienie pojawiają się w dystrybuancie rozkładu Gaussa – wybierane przy przeliczaniu między zmierzonym z-score a prawdopodobieństwem lub przy obliczaniu całki ogona rozkładu normalnego. Funkcje gamma i log-gamma pojawiają się w obliczeniach beta / chi-kwadrat / t-Studenta; gammaln jest formą stabilną numerycznie dla dużych argumentów, przy których sama gamma powodowałaby przepełnienie.

6.17.4. Liczby losowe

numpy.random dostarcza klasę Generator, która losuje próbki ze wspólnych rozkładów. Generator jest stanowy: każde wywołanie przesuwa jego wewnętrzny stan, więc kolejne wywołania zwracają niezależne próbki:

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

Wyjściowym typem danych jest zawsze float. size= przyjmuje liczbę całkowitą (wyjście jednowymiarowe) lub krotkę (wyjście n-wymiarowe); gdy pominięte, zwracana jest pojedyncza liczba float Pythona.

Generator nadaje się do symulacji, ditheringu, syntetycznych danych testowych i każdego innego zastosowania, w którym nie jest wymagana siła kryptograficzna. Nie nadaje się do kluczy ani tokenów; w tych przypadkach użyj systemowego źródła losowości przez os.

6.17.5. Dostępność w czasie kompilacji

To, czy każdy submoduł jest faktycznie obecny, zależy od sposobu, w jaki kamera została skompilowana. scipy.optimize i scipy.special nie są włączone na każdej kamerze; wywołanie funkcji, której kamera nie zawiera, zgłasza AttributeError. dir(sp), dir(sp.optimize), dir(np.random) i podobne raportują, co jest dostępne na docelowej kamerze.