6.17. Risolutori e numeri casuali

Quando la funzione in esame è definita da codice Python anziché da un buffer di campioni, la chiamata giusta è una famiglia di strumenti diversa: dov’è la radice della funzione, il suo minimo, il suo integrale su un dato intervallo? I sottomoduli scipy.integrate e scipy.optimize coprono questo lavoro. Ogni algoritmo richiama la funzione Python fornita dall’utente, quindi il costo per iterazione è più alto di una riduzione su buffer; la comodità sta nel non dover scrivere il risolutore.

Il sottomodulo scipy.special copre le funzioni speciali statistiche (funzione errore, gamma) che emergono quando si calcola la funzione di distribuzione cumulativa (CDF, la probabilità che un campione sia al più un dato valore) o la funzione di densità di probabilità (PDF, la verosimiglianza relativa a un dato valore) di una distribuzione di probabilità. numpy.random copre il generatore pseudo-casuale per dithering, simulazione e dati di test sintetici.

6.17.1. Integrazione numerica di una funzione richiamabile

Quando l’integrando è una funzione Python anziché un buffer di campioni, scipy.integrate espone quattro algoritmi di quadratura:

  • quad() – Gauss-Kronrod adattivo. L’impostazione predefinita giusta per integrandi regolari. Restituisce (value, error).

  • romberg() – Romberg / Newton-Cotes classico. Restituisce un singolo float. Deprecato a monte; incluso per compatibilità.

  • simpson() – regola di Simpson adattiva. Restituisce un singolo float.

  • tanhsinh() – quadratura doppio-esponenziale. Usala quando l’integrando ha singolarità agli estremi o un limite infinito. Restituisce (value, error).

L’integrale gaussiano valutato con la regola doppio-esponenziale (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))

Output:

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

6.17.2. Ricerca di radici e minimizzazione

scipy.optimize copre tre classici risolutori a variabile singola. Ogni iterazione richiama la funzione Python fornita dall’utente, quindi l’accelerazione rispetto a un risolutore in puro Python è modesta (circa 2x); la comodità sta nel non dover scrivere il risolutore.

  • bisect() – trova una radice di f su [a, b] dimezzando l’intervallo. f(a) e f(b) devono avere segni opposti:

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

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

  • newton() – trova una radice usando l’iterazione delle secanti / 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() – trova un minimo locale usando il metodo del simplesso discendente (Nelder-Mead):

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

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

Lo scope a variabile singola è sufficiente per la maggior parte delle ottimizzazioni lato camera – una costante di calibrazione di un sensore, il guadagno che massimizza una misura di contrasto, la soglia in cui la bimodalità di un istogramma è più netta. Per problemi a più variabili, la risposta giusta è di solito riformulare il problema come una piccola risoluzione di algebra lineare anziché ricorrere a un ottimizzatore non lineare generico.

6.17.3. Funzioni speciali

scipy.special espone una manciata di funzioni statistiche e di probabilità che si comportano come funzioni universali – accettano uno scalare, un iterabile o un ndarray e restituiscono un ndarray di float:

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 funzione errore e il suo complemento compaiono nella CDF di una gaussiana – l’applicazione di scelta per convertire tra uno z-score misurato e una probabilità o per calcolare l’integrale di coda di una distribuzione normale. Le funzioni gamma e log-gamma compaiono nei calcoli beta / chi-quadrato / student-t; gammaln è la forma numericamente stabile per argomenti grandi dove gamma stessa andrebbe in overflow.

6.17.4. Numeri casuali

numpy.random fornisce una classe Generator che estrae campioni da distribuzioni comuni. Il generatore è stateful: ogni chiamata fa avanzare il suo stato interno, quindi chiamate consecutive restituiscono campioni indipendenti:

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

Il dtype dell’output è sempre float. size= accetta un intero (output monodimensionale) o una tupla (output n-dimensionale); quando omesso, viene restituito un singolo float Python.

Il generatore è adatto a simulazione, dithering, dati di test sintetici e qualsiasi altra applicazione in cui non sia richiesta robustezza crittografica. Non è adatto a chiavi o token; per quelli usa la sorgente casuale di sistema attraverso os.

6.17.5. Disponibilità al momento della build

Se ciascun sottomodulo sia effettivamente presente dipende da come è stata costruita la camera. scipy.optimize e scipy.special non sono abilitati su ogni camera; chiamare una funzione che la camera non include solleva AttributeError. dir(sp), dir(sp.optimize), dir(np.random) e simili riportano ciò che è disponibile sulla camera di destinazione.