Tutorial SciPy — Calcolo Scientifico in Python

SciPy (Scientific Python) è una libreria open source che si basa su NumPy e aggiunge un'ampia raccolta di algoritmi per matematica, scienza e ingegneria. Mentre NumPy ti fornisce l'ndarray e le operazioni di base, SciPy ti offre i solutori specializzati: integrazione numerica, ottimizzazione, interpolazione, elaborazione dei segnali, statistica, algoritmi spaziali e altro ancora.

Questo capitolo tratta:

Installazione di SciPy e la convenzione di importazione
Algebra lineare (scipy.linalg)
Integrazione numerica e differenziazione (scipy.integrate)
Ottimizzazione — ricerca di minimi e radici (scipy.optimize)
Interpolazione (scipy.interpolate)
Statistica (scipy.stats)
Elaborazione di immagini N-dimensionali (scipy.ndimage)

Installazione di SciPy

SciPy è incluso nella distribuzione Anaconda. Per installarlo con pip:

pip install scipy

SciPy dipende da NumPy, che pip installa automaticamente se non è già presente.

Convenzione di importazione

SciPy è organizzato in sotto-pacchetti. Importa solo i sotto-pacchetti di cui hai bisogno anziché l'intera libreria:

import numpy as np
from scipy import linalg, integrate, optimize, interpolate, stats, ndimage

Puoi anche verificare la versione installata:

import scipy
print(scipy.__version__)  # e.g. 1.13.0

Algebra Lineare

scipy.linalg estende le routine di algebra lineare di NumPy con decomposizioni e solutori aggiuntivi. Tutte le funzioni operano su array NumPy standard.

Determinante e inversa

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

det = linalg.det(a)
print(det)          # -2.0

inv = linalg.inv(a)
print(inv)
# [[-2.   1. ]
#  [ 1.5 -0.5]]

det([[1,2],[3,4]]) = 1×4 − 2×3 = −2. L'inversa soddisfa a @ inv == I.

Autovalori e autovettori

Gli autovalori descrivono come una matrice deforma lo spazio; gli autovettori indicano le direzioni che non ruotano.

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

eigenvalues, eigenvectors = linalg.eig(a)
print(eigenvalues)
# [-0.37228132+0.j  5.37228132+0.j]
print(eigenvectors)
# [[-0.82456484 -0.41597356]
#  [ 0.56576746 -0.90937671]]

Ogni colonna di eigenvectors corrisponde all'autovalore corrispondente. Gli autovalori vengono restituiti come numeri complessi anche quando la parte immaginaria è zero.

Decomposizione ai Valori Singolari (SVD)

La SVD fattorizza una matrice A in tre matrici U, s, Vt tali che A = U @ diag(s) @ Vt. È alla base dell'analisi delle componenti principali (PCA) e di molte tecniche di riduzione della dimensionalità.

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

u, s, vt = linalg.svd(a)
print(u)
# [[-0.40455358 -0.9145143 ]
#  [-0.9145143   0.40455358]]
print(s)    # [5.4649857  0.36596619]
print(vt)
# [[-0.57604844 -0.81741556]
#  [ 0.81741556 -0.57604844]]

Soluzione di un sistema lineare

linalg.solve è il metodo corretto per risolvere Ax = b. È più veloce e numericamente più stabile rispetto al calcolo dell'inversa seguito da una moltiplicazione.

import numpy as np
from scipy import linalg

# Solve: 1x + 2y = 5
#        3x + 4y = 11
A = np.array([[1, 2],
              [3, 4]])
b = np.array([5, 11])

x = linalg.solve(A, b)
print(x)   # [1. 2.]

# Verify: A @ x should equal b
print(np.allclose(A @ x, b))  # True

Integrazione Numerica

scipy.integrate fornisce routine per il calcolo di integrali definiti quando una soluzione analitica è impraticabile.

Integrazione a variabile singola con `quad`

integrate.quad utilizza la quadratura adattiva per integrare una funzione su un intervallo. Restituisce il risultato e una stima dell'errore assoluto.

import numpy as np
from scipy import integrate

# Integrate f(x) = x^2 + 2x + 1 from 0 to 1
# Analytical result: [x^3/3 + x^2 + x] from 0 to 1 = 1/3 + 1 + 1 = 7/3
def f(x):
    return x**2 + 2*x + 1

result, error = integrate.quad(f, 0, 1)
print(result)   # 2.3333333333333335
print(error)    # ~2.6e-14  (absolute error estimate)

Differenziazione numerica

approx_fprime di SciPy calcola un gradiente alle differenze finite. Per funzioni scalari, la derivative alle differenze centrali da scipy.misc è più semplice:

import numpy as np
from scipy.optimize import approx_fprime

# Derivative of sin(x) at x = 0 should be cos(0) = 1
result = approx_fprime([0.0], lambda x: np.sin(x[0]), 1e-8)
print(result[0])   # ~1.0

Ottimizzazione

scipy.optimize trova minimi, massimi (minimizzando il negativo) e radici delle funzioni.

Minimizzazione di una funzione multivariata

optimize.minimize supporta molti metodi (Nelder-Mead, BFGS, L-BFGS-B, …). Il metodo predefinito viene scelto automaticamente.

from scipy import optimize

# Minimize f(x) = x^2 + 2x + 1 = (x + 1)^2  — minimum at x = -1
def f(x):
    return x**2 + 2*x + 1

result = optimize.minimize(f, x0=0)   # x0 is the starting guess
print(result.success)   # True
print(result.x)         # [-1.00000001]  (near -1)
print(result.fun)       # ~0.0           (minimum value)

Minimizzazione di una funzione scalare su un intervallo limitato

optimize.minimize_scalar è più semplice per i problemi a variabile singola. Fornisci sempre bounds con method='bounded' quando la funzione non ha un minimo globale (cioè è illimitata):

from scipy import optimize

# Minimize h(x) = x^2 - 4x + 3 over [0, 4]  — minimum at x = 2, h(2) = -1
def h(x):
    return x**2 - 4*x + 3

result = optimize.minimize_scalar(h, bounds=(0, 4), method='bounded')
print(result.x)     # ~2.0
print(result.fun)   # -1.0

Ricerca delle radici

optimize.root_scalar trova dove una funzione attraversa lo zero:

from scipy.optimize import root_scalar

# Solve x^2 - 4 = 0 in the interval [0, 3]  — root at x = 2
res = root_scalar(lambda x: x**2 - 4, bracket=[0, 3])
print(res.root)     # 2.0

Interpolazione

scipy.interpolate adatta una curva regolare a punti dati in modo da poter stimare i valori intermedi.

Interpolazione 1-D

interp1d crea una funzione di interpolazione richiamabile a partire da coppie discrete (x, y). Il parametro kind seleziona il metodo: 'linear' (predefinito), 'quadratic' o 'cubic'.

import numpy as np
from scipy.interpolate import interp1d

# Sample points from y = x^2
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

f_linear = interp1d(x, y)           # piecewise linear
f_cubic  = interp1d(x, y, kind='cubic')  # cubic spline

# Estimate y at x = 2.5  (exact value: 2.5^2 = 6.25)
print(float(f_linear(2.5)))   # 6.5   (linear — slightly off)
print(float(f_cubic(2.5)))    # 6.25  (cubic — matches exactly for polynomials)

L'interpolazione cubica recupera il risultato esatto perché i dati provengono da un polinomio quadratico e la spline cubica è abbastanza flessibile da adattarsi perfettamente ad esso.

Statistica

scipy.stats contiene oltre 80 distribuzioni di probabilità continue e discrete, più una raccolta di test statistici.

Statistiche descrittive e distribuzione normale

import numpy as np
from scipy.stats import norm

# Standard normal distribution (mean=0, std=1)
print(norm.pdf(0))     # 0.3989422804014327  — probability density at x = 0
print(norm.cdf(1.96))  # 0.9750021048517795  — P(X <= 1.96)
print(norm.ppf(0.975)) # 1.9599639845400536  — inverse CDF (quantile function)

# Fit a normal distribution to data
data = np.array([2.1, 3.3, 2.8, 3.1, 2.5, 3.0, 2.7])
mu, sigma = norm.fit(data)
print(f'Fitted mean: {mu:.4f}, std: {sigma:.4f}')

Test di ipotesi

scipy.stats include t-test, test chi-quadrato, ANOVA, test di Kolmogorov-Smirnov e altri ancora.

import numpy as np
from scipy.stats import ttest_1samp, ttest_ind

# One-sample t-test: is the sample mean significantly different from 5?
sample = np.array([4.8, 5.1, 4.9, 5.3, 5.2, 4.7, 5.0])
t_stat, p_value = ttest_1samp(sample, popmean=5.0)
print(f't = {t_stat:.4f}, p = {p_value:.4f}')
# p > 0.05 → no significant difference from 5

# Two-sample t-test: are these two groups different?
group_a = np.array([5.1, 5.3, 4.9, 5.2, 5.0])
group_b = np.array([6.1, 6.3, 5.8, 6.0, 5.9])
t_stat2, p_value2 = ttest_ind(group_a, group_b)
print(f't = {t_stat2:.4f}, p = {p_value2:.6f}')
# Very small p → groups are significantly different

Variabili casuali e campionamento

Ogni distribuzione in scipy.stats espone la stessa interfaccia: pdf, cdf, ppf, rvs (variabili casuali) e fit.

from scipy.stats import norm, poisson

# Draw 5 samples from a normal distribution with mean=10, std=2
samples = norm.rvs(loc=10, scale=2, size=5, random_state=42)
print(samples.round(2))   # [10.99  9.72 11.3  13.05  9.53]

# Poisson distribution: P(X = k) for mean lambda=3
for k in range(6):
    print(f'P(X={k}) = {poisson.pmf(k, mu=3):.4f}')

Elaborazione di Immagini N-Dimensionali

scipy.ndimage opera su array di qualsiasi dimensionalità (immagini, volumi, cubi di serie temporali). Di seguito è riportato un esempio autonomo che utilizza un array 2-D sintetico, senza necessità di file immagine esterni.

Sfocatura gaussiana ed etichettatura delle regioni connesse

import numpy as np
from scipy import ndimage

# Create a synthetic 5x5 "image" with a bright spot in the centre
image = np.array([
    [0, 0, 0, 0, 0],
    [0, 1, 1, 1, 0],
    [0, 1, 5, 1, 0],
    [0, 1, 1, 1, 0],
    [0, 0, 0, 0, 0],
], dtype=float)

# Smooth the image with a Gaussian filter (sigma controls the blur radius)
blurred = ndimage.gaussian_filter(image, sigma=1)
print('Center value after blur:', round(blurred[2, 2], 4))
# 1.4166  — the bright peak is spread across neighbouring pixels

# Label connected non-zero regions (like counting distinct objects)
binary = image > 0
labeled, num_regions = ndimage.label(binary)
print('Number of connected regions:', num_regions)  # 1
print(labeled)
# [[0 0 0 0 0]
#  [0 1 1 1 0]
#  [0 1 1 1 0]
#  [0 1 1 1 0]
#  [0 0 0 0 0]]

Operazioni comuni di `ndimage`

Funzione	Cosa fa
`gaussian_filter(a, sigma)`	Applica un filtro gaussiano per smussare
`sobel(a)`	Rileva i bordi (gradiente di Sobel)
`label(a)`	Etichetta le regioni connesse in un array binario
`binary_dilation(a)`	Espande le regioni in primo piano
`zoom(a, factor)`	Ridimensiona un array
`rotate(a, angle)`	Ruota un array (in gradi)

Quando Usare SciPy Rispetto ad Altre Librerie

Attività	Strumento preferito
Creazione di array e operazioni matematiche di base	NumPy
DataFrame, serie temporali, I/O	Pandas
Algoritmi scientifici (integrazione, ottimizzazione)	SciPy
Machine learning	scikit-learn (basato su SciPy)
Visualizzazione	Matplotlib

SciPy non è un sostituto di NumPy — dipende da NumPy e lo estende. In pratica importerai entrambi.

Riferimento Rapido

Sotto-pacchetto	Funzioni principali
`scipy.linalg`	`det`, `inv`, `eig`, `svd`, `solve`
`scipy.integrate`	`quad`, `dblquad`, `solve_ivp`
`scipy.optimize`	`minimize`, `minimize_scalar`, `root_scalar`
`scipy.interpolate`	`interp1d`, `CubicSpline`, `griddata`
`scipy.stats`	`norm`, `ttest_1samp`, `ttest_ind`, `chi2_contingency`
`scipy.ndimage`	`gaussian_filter`, `label`, `sobel`, `zoom`
`scipy.signal`	`butter`, `lfilter`, `find_peaks`
`scipy.spatial`	`distance`, `KDTree`, `ConvexHull`