SciPy: Cálculo científico avanzado en Python

SciPy es la biblioteca de cálculo científico que amplía NumPy con algoritmos numéricos contrastados: optimización, ajuste de curvas, estadística, integración de ecuaciones diferenciales, interpolación, procesamiento de señales y matrices dispersas. Donde NumPy proporciona la estructura de datos (el ndarray) y las operaciones elementales, SciPy proporciona los métodos numéricos de alto nivel, en su mayoría respaldados por bibliotecas de referencia en Fortran y C (LAPACK, MINPACK, FFTPACK). Esta guía recorre los submódulos más utilizados con ejemplos completos y ejecutables.

SciPy frente a NumPy

La frontera es fácil de recordar:

• NumPy: el array ndarray, la aritmética vectorizada, el álgebra lineal básica, la generación de números aleatorios.
• SciPy: los algoritmos que operan sobre esos arrays, con implementaciones robustas fruto de décadas de cálculo científico.

SciPy depende de NumPy y se importa submódulo a submódulo:

import numpy as np
from scipy import optimize, stats, integrate, interpolate, signal, sparse

# Verificar la versión instalada
import scipy
print(scipy.__version__)     # 1.14.1 (por ejemplo)

Existen algunos solapamientos (scipy.linalg amplía numpy.linalg con descomposiciones adicionales; scipy.fft sustituye con ventaja a numpy.fft), y en esos casos la versión de SciPy suele ser más completa.

scipy.optimize: minimización y ajuste de curvas

Minimizar una función con minimize

minimize encuentra un mínimo local de una función escalar de varias variables.

from scipy.optimize import minimize

# Función de Rosenbrock, un banco de pruebas clásico
def rosenbrock(p):
    x, y = p
    return (1 - x) ** 2 + 100 * (y - x ** 2) ** 2

resultat = minimize(rosenbrock, x0=[0.0, 0.0], method="Nelder-Mead")

print(resultat.success)      # True
print(resultat.x)            # [0.99998736 0.99997369] -> minimum en (1, 1)
print(resultat.fun)          # 5.9e-10 (valeur au minimum, quasi nulle)
print(resultat.nit)          # nombre d'itérations

La elección del método depende del problema: "Nelder-Mead" no requiere gradiente, "BFGS" (por defecto sin restricciones) lo estima mediante diferencias finitas, "L-BFGS-B" admite límites, "SLSQP" gestiona restricciones generales.

# Minimisation avec bornes : x dans [0, 2], y dans [0, 2]
res = minimize(rosenbrock, x0=[0.5, 0.5], method="L-BFGS-B",
               bounds=[(0, 2), (0, 2)])
print(res.x)                 # [1. 1.]

Para encontrar la raíz de una ecuación en lugar de un mínimo, utilice optimize.root_scalar o optimize.brentq.

Ajustar un modelo a los datos con curve_fit

curve_fit estima los parámetros de un modelo mediante mínimos cuadrados no lineales. Ejemplo: ajustar una decaimiento exponencial a mediciones con ruido.

from scipy.optimize import curve_fit

# Modèle : la première variable est toujours x, les suivantes sont les paramètres
def modele(x, a, tau, c):
    return a * np.exp(-x / tau) + c

# Données synthétiques bruitées
rng = np.random.default_rng(42)
x_donnees = np.linspace(0, 10, 50)
y_vrai = modele(x_donnees, a=3.0, tau=2.5, c=0.5)
y_donnees = y_vrai + rng.normal(0, 0.15, size=x_donnees.size)

# Ajustement (p0 : estimation initiale des paramètres)
params, covariance = curve_fit(modele, x_donnees, y_donnees, p0=[1, 1, 0])

a_est, tau_est, c_est = params
incertitudes = np.sqrt(np.diag(covariance))

print(f"a   = {a_est:.3f} ± {incertitudes[0]:.3f}")   # a   = 3.011 ± 0.052
print(f"tau = {tau_est:.3f} ± {incertitudes[1]:.3f}") # tau = 2.487 ± 0.098
print(f"c   = {c_est:.3f} ± {incertitudes[2]:.3f}")   # c   = 0.503 ± 0.031

La diagonal de la matriz de covarianza da la varianza de cada parámetro estimado: su raíz cuadrada proporciona una incertidumbre de una desviación típica. Proporcione siempre un p0 razonable en los modelos no lineales, a riesgo de que converjan hacia un óptimo equivocado.

scipy.stats: distribuciones y pruebas estadísticas

Trabajar con distribuciones

Cada distribución expone la misma interfaz: pdf (densidad), cdf (función de distribución), ppf (cuantiles, inversa de la cdf) y rvs (muestreo).

from scipy import stats

# Loi normale de moyenne 100 et d'écart-type 15
loi = stats.norm(loc=100, scale=15)

print(loi.pdf(100))          # 0.02659... (densité au sommet)
print(loi.cdf(115))          # 0.8413... (P(X <= 115))
print(loi.ppf(0.975))        # 129.399... (quantile à 97,5 %)
print(loi.rvs(size=3, random_state=42))   # [107.45  84.44 109.72]

# Autres lois du même moule : stats.t, stats.chi2, stats.binom,
# stats.poisson, stats.expon, stats.uniform...

Prueba t: comparar dos medias

rng = np.random.default_rng(1)
groupe_a = rng.normal(loc=52, scale=8, size=40)   # ex. scores avec méthode A
groupe_b = rng.normal(loc=48, scale=8, size=40)   # ex. scores avec méthode B

t_stat, p_valeur = stats.ttest_ind(groupe_a, groupe_b)
print(f"t = {t_stat:.3f}, p = {p_valeur:.4f}")
# t = 2.279, p = 0.0255 (par exemple)

if p_valeur < 0.05:
    print("Différence significative au seuil de 5 %")
# Variante appariée (mesures avant/après sur les mêmes sujets) : stats.ttest_rel
# Variances inégales : ttest_ind(..., equal_var=False) (test de Welch)

Prueba de chi-cuadrado: independencia de dos variables categóricas

# Tableau de contingence : lignes = formation suivie, colonnes = réussite
observations = np.array([[45, 15],    # avec formation : 45 réussites, 15 échecs
                         [30, 30]])   # sans formation : 30 réussites, 30 échecs

chi2, p, ddl, attendus = stats.chi2_contingency(observations)
print(f"chi2 = {chi2:.3f}, p = {p:.4f}, ddl = {ddl}")
# chi2 = 6.806, p = 0.0091, ddl = 1
print(attendus)              # effectifs attendus sous l'hypothèse d'indépendance
# [[37.5 22.5]
#  [37.5 22.5]]

Correlaciones

x = rng.normal(size=100)
y = 0.7 * x + rng.normal(scale=0.5, size=100)

r, p = stats.pearsonr(x, y)          # corrélation linéaire
print(f"Pearson  r = {r:.3f}, p = {p:.2e}")   # r = 0.812, p = 1.5e-24

rho, p = stats.spearmanr(x, y)       # corrélation de rangs (relations monotones)
print(f"Spearman rho = {rho:.3f}")   # rho = 0.795

Utilice Spearman cuando la relación sea monótona pero no lineal, o en presencia de valores atípicos.

scipy.integrate: integrales y ecuaciones diferenciales

Integración numérica con quad

from scipy.integrate import quad

# Intégrale de exp(-x²) entre 0 et l'infini (vaut sqrt(pi)/2)
valeur, erreur = quad(lambda x: np.exp(-x ** 2), 0, np.inf)
print(valeur)                        # 0.8862269254527579
print(np.sqrt(np.pi) / 2)            # 0.8862269254527579
print(erreur)                        # 7.1e-09 (estimation de l'erreur)

Resolver una EDO con solve_ivp

Ejemplo: el oscilador armónico amortiguado, descrito por la ecuación de segundo orden x” + 2ζω x’ + ω² x = 0, reescrita como un sistema de primer orden.

from scipy.integrate import solve_ivp

omega = 2.0    # pulsation propre
zeta = 0.15    # taux d'amortissement

def oscillateur(t, y):
    x, v = y                              # position, vitesse
    return [v, -2 * zeta * omega * v - omega ** 2 * x]

solution = solve_ivp(
    oscillateur,
    t_span=(0, 20),                       # intervalle d'intégration
    y0=[1.0, 0.0],                        # position 1, vitesse nulle à t=0
    t_eval=np.linspace(0, 20, 400),       # points où évaluer la solution
    method="RK45",                        # Runge-Kutta adaptatif (défaut)
)

print(solution.success)                   # True
print(solution.y.shape)                   # (2, 400) -> position et vitesse
print(solution.y[0, :3])                  # [1.  0.995 0.980] (position amortie)

# Tracé de la position au cours du temps
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(solution.t, solution.y[0])
ax.set_xlabel("Temps (s)")
ax.set_ylabel("Position")
ax.set_title("Oscillateur harmonique amorti")
plt.show()

Para los sistemas rígidos (constantes de tiempo muy distintas), utilice method="Radau" o method="BDF".

scipy.interpolate: interpolación de datos

La interpolación construye una función continua que pasa por puntos de medición discretos.

from scipy.interpolate import CubicSpline, interp1d

# Points de mesure épars
x_mesures = np.array([0, 1, 2, 3, 4, 5], dtype=float)
y_mesures = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

# Spline cubique : lisse, dérivable, recommandée
spline = CubicSpline(x_mesures, y_mesures)

x_fin = np.linspace(0, 5, 100)
y_interpole = spline(x_fin)

print(spline(2.5))                   # 0.5546... (valeur interpolée entre 2 et 3)
print(spline(2.5, 1))                # dérivée première au même point

# Interpolation linéaire simple
lineaire = interp1d(x_mesures, y_mesures, kind="linear")
print(lineaire(2.5))                 # 0.5 (moyenne des points voisins)

Atención: la interpolación solo es fiable dentro del dominio de las mediciones. Para extrapolar o para suavizar datos con ruido, es preferible un ajuste de modelo (curve_fit) o una spline de suavizado (scipy.interpolate.make_smoothing_spline).

scipy.signal y scipy.fft: procesamiento de señales

Análisis espectral con la FFT

from scipy.fft import fft, fftfreq

# Signal échantillonné : somme de deux sinusoïdes (5 Hz et 50 Hz) plus du bruit
fe = 500                              # fréquence d'échantillonnage (Hz)
duree = 2.0
t = np.linspace(0, duree, int(fe * duree), endpoint=False)
rng = np.random.default_rng(7)
signal_brut = (np.sin(2 * np.pi * 5 * t)
               + 0.5 * np.sin(2 * np.pi * 50 * t)
               + 0.3 * rng.normal(size=t.size))

# Transformée de Fourier
spectre = fft(signal_brut)
frequences = fftfreq(t.size, d=1 / fe)

# Amplitudes sur les fréquences positives
moitie = t.size // 2
amplitudes = 2 / t.size * np.abs(spectre[:moitie])

pics = frequences[:moitie][amplitudes > 0.3]
print(pics)                           # [ 5. 50.] -> les deux fréquences retrouvées

Filtrado con un filtro de Butterworth

from scipy.signal import butter, filtfilt

# Filtre passe-bas d'ordre 4, fréquence de coupure 20 Hz
b, a = butter(N=4, Wn=20, btype="low", fs=fe)

# filtfilt applique le filtre dans les deux sens : aucun déphasage
signal_filtre = filtfilt(b, a, signal_brut)

# La composante à 50 Hz et le bruit haute fréquence sont supprimés,
# la sinusoïde à 5 Hz est préservée.

scipy.signal ofrece también la detección de picos (find_peaks), los espectrogramas (ShortTimeFFT), la convolución y el diseño de numerosos tipos de filtros (Chebyshev, elípticos, FIR).

from scipy.signal import find_peaks

pics_indices, proprietes = find_peaks(signal_filtre, height=0.5, distance=fe // 10)
print(len(pics_indices))              # 10 -> un pic par période de la sinusoïde 5 Hz

scipy.sparse: matrices dispersas

Cuando una matriz contiene una abrumadora mayoría de ceros (grafos, sistemas de ecuaciones procedentes de mallados, matrices término-documento), almacenarla de forma densa desperdicia memoria y cálculo.

from scipy import sparse

# Matrice creuse au format COO (construction) puis CSR (calcul)
lignes = np.array([0, 1, 2, 2])
colonnes = np.array([1, 2, 0, 2])
valeurs = np.array([3.0, 4.0, 5.0, 6.0])

m = sparse.coo_array((valeurs, (lignes, colonnes)), shape=(3, 3)).tocsr()
print(m.toarray())
# [[0. 3. 0.]
#  [0. 0. 4.]
#  [5. 0. 6.]]

print(m.nnz)                          # 4 éléments non nuls stockés (sur 9)

# Produit matrice-vecteur : seuls les éléments non nuls sont parcourus
v = np.array([1.0, 2.0, 3.0])
print(m @ v)                          # [ 6. 12. 23.]

# Résolution d'un système linéaire creux
from scipy.sparse.linalg import spsolve
b = np.array([6.0, 12.0, 23.0])
x = spsolve(m.tocsc(), b)
print(x)                              # [1. 2. 3.]

Órdenes de magnitud: una matriz densa de 100 000 × 100 000 en float64 ocuparía 80 GB; la misma matriz con un 0,01 % de elementos no nulos cabe en unos pocos megabytes en formato CSR.

Tabla resumen de los submódulos

Submódulo	Ámbito	Funciones destacadas
scipy.optimize	Optimización, ajuste, raíces	minimize, curve_fit, root_scalar, linprog
scipy.stats	Distribuciones, pruebas, correlaciones	norm, ttest_ind, chi2_contingency, pearsonr
scipy.integrate	Integrales, ecuaciones diferenciales	quad, solve_ivp
scipy.interpolate	Interpolación, splines	CubicSpline, interp1d, make_smoothing_spline
scipy.signal	Filtrado, picos, espectrogramas	butter, filtfilt, find_peaks
scipy.fft	Transformadas de Fourier	fft, ifft, fftfreq, rfft
scipy.sparse	Matrices dispersas	csr_array, coo_array, sparse.linalg.spsolve
scipy.linalg	Álgebra lineal ampliada	lu, svd, expm, solve_banded
scipy.spatial	Geometría, distancias	KDTree, distance.cdist, ConvexHull
scipy.ndimage	Procesamiento de imágenes n-D	gaussian_filter, label, rotate

Buenas prácticas

• Importe los submódulos explícitamente (from scipy import optimize): import scipy por sí solo no carga los submódulos.
• Compruebe siempre resultat.success después de minimize o solve_ivp: un algoritmo que no ha convergido devuelve igualmente un resultado.
• Proporcione una estimación inicial p0 razonable a curve_fit: los mínimos cuadrados no lineales convergen hacia un óptimo local que depende del punto de partida.
• Interprete los valores p con prudencia: un valor p inferior a 0,05 no mide ni el tamaño del efecto ni su importancia práctica; indique también las medias y los intervalos de confianza.
• Elija el método adecuado al problema: Radau o BDF para las EDO rígidas, L-BFGS-B en cuanto haya límites, la prueba de Welch cuando las varianzas difieren.
• No extrapole con una interpolación: fuera del dominio de las mediciones, una spline cúbica diverge rápidamente.
• Pase a las matrices dispersas en cuanto la densidad de elementos no nulos baje de unos pocos puntos porcentuales en matrices grandes.

Conclusión

SciPy convierte Python en un entorno de cálculo científico completo: optimización y ajuste de modelos con optimize, inferencia estadística con stats, resolución de ecuaciones diferenciales con integrate, reconstrucción de funciones con interpolate, análisis en frecuencia y filtrado con signal y fft, y escalado con sparse. Cada submódulo encapsula algoritmos de referencia validados desde hace décadas, lo que le evita reimplementar métodos numéricos delicados. Combinado con NumPy para las estructuras de datos, Pandas para las tablas etiquetadas y Matplotlib para la visualización, SciPy completa la caja de herramientas estándar de cualquier trabajo científico o de ingeniería en Python.