6.17. Resolvedores y números aleatorios

Cuando la función bajo estudio está definida por código Python en lugar de por un búfer de muestras, la opción adecuada es una familia de herramientas diferente: ¿dónde está la raíz de la función, su mínimo, su integral sobre un intervalo dado? Los submódulos scipy.integrate y scipy.optimize cubren ese trabajo. Cada algoritmo vuelve a llamar a la función Python suministrada por el usuario, por lo que el coste por iteración es mayor que el de una reducción sobre un búfer; la comodidad está en no tener que escribir el resolvedor.

El submódulo scipy.special cubre las funciones especiales estadísticas (función de error, gamma) que surgen al calcular la función de distribución acumulada (CDF, la probabilidad de que una muestra sea como máximo un valor dado) o la función de densidad de probabilidad (PDF, la verosimilitud relativa en un valor dado) de una distribución de probabilidad. numpy.random cubre el generador pseudoaleatorio para difuminado, simulación y datos de prueba sintéticos.

6.17.1. Integración numérica de un invocable

Cuando el integrando es una función Python en lugar de un búfer de muestras, scipy.integrate expone cuatro algoritmos de cuadratura:

  • quad() – Gauss-Kronrod adaptativa. La opción por defecto adecuada para integrandos suaves. Devuelve (value, error).

  • romberg() – Romberg / Newton-Cotes clásica. Devuelve un único float. Obsoleta en el proyecto original; incluida por compatibilidad.

  • simpson() – regla de Simpson adaptativa. Devuelve un único float.

  • tanhsinh() – cuadratura doblemente exponencial. Úsela cuando el integrando tenga singularidades en los extremos o un límite infinito. Devuelve (value, error).

La integral gaussiana evaluada con la regla doblemente exponencial (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))

Salida:

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

6.17.2. Búsqueda de raíces y minimización

scipy.optimize cubre tres resolvedores clásicos de una sola variable. Cada iteración vuelve a llamar a la función Python suministrada por el usuario, por lo que la aceleración respecto a un resolvedor en Python puro es modesta (aproximadamente 2x); la comodidad está en no tener que escribir el resolvedor.

  • bisect() – encuentra una raíz de f en [a, b] reduciendo el intervalo a la mitad. f(a) y f(b) deben tener signos opuestos:

    def f(x):
        return x * x - 1
    
    sp.optimize.bisect(f, 0, 4)        # ~1.0
    
  • newton() – encuentra una raíz mediante iteración por secante / Newton-Raphson:

    def f(x):
        return x * x * x - 2.0
    
    sp.optimize.newton(f, 3., tol=0.001, rtol=0.01)
    # ~1.260
    
  • fmin() – encuentra un mínimo local usando el método del símplex descendente (Nelder-Mead):

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

El alcance de una sola variable basta para la mayoría de las optimizaciones del lado de la cámara – una constante de calibración de un sensor, la ganancia que maximiza una medida de contraste, el umbral en el que la bimodalidad de un histograma es más nítida. Para problemas de varias variables, la respuesta correcta suele ser reformular el problema como una pequeña resolución de álgebra lineal en lugar de recurrir a un optimizador no lineal general.

6.17.3. Funciones especiales

scipy.special expone un puñado de funciones estadísticas y de probabilidad que se comportan como funciones universales – aceptan un escalar, un iterable o un ndarray y devuelven un ndarray de coma flotante:

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

La función de error y su complemento aparecen en la CDF de una gaussiana – la aplicación idónea para convertir entre una puntuación z medida y una probabilidad, o para calcular la integral de cola de una distribución normal. Las funciones gamma y log-gamma aparecen en cálculos de beta / chi cuadrado / t de Student; gammaln es la forma numéricamente estable para argumentos grandes en los que el propio gamma desbordaría.

6.17.4. Números aleatorios

numpy.random proporciona una clase Generator que extrae muestras de distribuciones comunes. El generador tiene estado: cada llamada avanza su estado interno, de modo que llamadas consecutivas devuelven muestras independientes:

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

El dtype de salida es siempre float. size= acepta un entero (salida unidimensional) o una tupla (salida n-dimensional); cuando se omite, se devuelve un único float de Python.

El generador es adecuado para simulación, difuminado, datos de prueba sintéticos y cualquier otra aplicación donde no se requiera fortaleza criptográfica. No es adecuado para claves ni tokens; para eso, use la fuente aleatoria del sistema a través de os.

6.17.5. Disponibilidad en tiempo de compilación

Que cada submódulo esté realmente presente depende de cómo se compiló la cámara. scipy.optimize y scipy.special no están habilitados en todas las cámaras; llamar a una función que la cámara no incluye lanza AttributeError. dir(sp), dir(sp.optimize), dir(np.random) y similares informan de lo que está disponible en la cámara objetivo.