6.17. Решатели и случайные числа

Когда исследуемая функция задана кодом на Python, а не буфером отсчётов, на помощь приходит другое семейство инструментов: где находится корень функции, её минимум, её интеграл на заданном интервале? Эту работу покрывают подмодули scipy.integrate и scipy.optimize. Каждый алгоритм обращается обратно к предоставленной пользователем функции Python, поэтому стоимость одной итерации выше, чем у редукции буфера; удобство в том, что не приходится писать сам решатель.

Подмодуль scipy.special покрывает статистические специальные функции (функция ошибок, гамма-функция), которые возникают при вычислении функции распределения (CDF, вероятность того, что отсчёт не превышает заданного значения) или функции плотности вероятности (PDF, относительная вероятность при заданном значении) для распределения вероятностей. numpy.random покрывает генератор псевдослучайных чисел для дизеринга, моделирования и синтетических тестовых данных.

6.17.1. Численное интегрирование вызываемого объекта

Когда подынтегральная функция – это функция Python, а не буфер отсчётов, scipy.integrate предоставляет четыре алгоритма квадратуры:

  • quad() – адаптивный метод Гаусса-Кронрода. Правильный выбор по умолчанию для гладких подынтегральных функций. Возвращает (value, error).

  • romberg() – классический метод Ромберга / Ньютона-Котеса. Возвращает одно число с плавающей точкой. Объявлен устаревшим в апстриме; включён для совместимости.

  • simpson() – адаптивный метод Симпсона. Возвращает одно число с плавающей точкой.

  • 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 покрывает три классических одномерных решателя. Каждая итерация обращается обратно к предоставленной пользователем функции Python, поэтому ускорение по сравнению с чисто Python-решателем умеренное (примерно в 2 раза); удобство в том, что не приходится писать сам решатель.

  • bisect() – находит корень f на [a, b] путём деления интервала пополам. 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 и возвращают 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-показателем и вероятностью или для вычисления хвостового интеграла нормального распределения. Гамма-функция и логарифм гамма-функции возникают в вычислениях бета-, хи-квадрат- и стьюдентова распределений; 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))

Тип данных вывода всегда float. size= принимает целое число (одномерный вывод) или кортеж (n-мерный вывод); если опущен, возвращается одно число Python с плавающей точкой.

Генератор подходит для моделирования, дизеринга, синтетических тестовых данных и любого другого приложения, где не требуется криптографическая стойкость. Он не подходит для ключей или токенов; для них используйте системный источник случайности через os.

6.17.5. Доступность во время сборки

Присутствует ли каждый подмодуль на самом деле, зависит от того, как была собрана камера. scipy.optimize и scipy.special включены не на каждой камере; вызов функции, которой камера не содержит, вызывает исключение AttributeError. dir(sp), dir(sp.optimize), dir(np.random) и подобные сообщают, что доступно на целевой камере.