Tutorial de SciPy — Computación Científica en Python

SciPy (Scientific Python) es una biblioteca de código abierto que se basa en NumPy y añade una gran colección de algoritmos para matemáticas, ciencias e ingeniería. Mientras NumPy te proporciona el ndarray y las operaciones básicas, SciPy te ofrece los solucionadores especializados: integración numérica, optimización, interpolación, procesamiento de señales, estadística, algoritmos espaciales y más.

Este capítulo cubre:

Instalación de SciPy y la convención de importación
Álgebra lineal (scipy.linalg)
Integración y diferenciación numéricas (scipy.integrate)
Optimización — búsqueda de mínimos y raíces (scipy.optimize)
Interpolación (scipy.interpolate)
Estadística (scipy.stats)
Procesamiento de imágenes N-dimensional (scipy.ndimage)

Instalación de SciPy

SciPy está incluido en la distribución Anaconda. Para instalarlo con pip:

pip install scipy

SciPy depende de NumPy, que pip instala automáticamente si aún no está presente.

Convención de importación

SciPy está organizado en sub-paquetes. Importa únicamente los sub-paquetes que necesites en lugar de toda la biblioteca:

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

También puedes comprobar la versión instalada:

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

Álgebra Lineal

scipy.linalg extiende las rutinas de álgebra lineal de NumPy con descomposiciones y solucionadores adicionales. Todas las funciones operan sobre arrays estándar de NumPy.

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. La inversa satisface a @ inv == I.

Valores propios y vectores propios

Los valores propios describen cómo una matriz estira el espacio; los vectores propios indican las direcciones que no se rotan.

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

Cada columna de eigenvectors corresponde al valor propio correspondiente. Los valores propios se devuelven como números complejos incluso cuando la parte imaginaria es cero.

Descomposición en Valores Singulares (SVD)

SVD factoriza una matriz A en tres matrices U, s, Vt de forma que A = U @ diag(s) @ Vt. Es la base del análisis de componentes principales (PCA) y muchas técnicas de reducción de dimensionalidad.

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

Resolución de un sistema lineal

linalg.solve es la forma correcta de resolver Ax = b. Es más rápida y numéricamente más estable que calcular la inversa y multiplicar.

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

Integración Numérica

scipy.integrate proporciona rutinas para calcular integrales definidas cuando una solución analítica no es práctica.

Integración de una variable con `quad`

integrate.quad utiliza cuadratura adaptativa para integrar una función sobre un intervalo. Devuelve el resultado y una estimación del error absoluto.

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)

Diferenciación numérica

La función approx_fprime de SciPy calcula un gradiente por diferencias finitas. Para funciones escalares, la derivative de diferencias centrales de scipy.misc es más sencilla:

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

Optimización

scipy.optimize encuentra mínimos, máximos (minimizando el negativo) y raíces de funciones.

Minimizar una función multivariante

optimize.minimize admite muchos métodos (Nelder-Mead, BFGS, L-BFGS-B, …). El método predeterminado se elige automáticamente.

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)

Minimizar una función escalar en un intervalo acotado

optimize.minimize_scalar es más sencillo para problemas de una sola variable. Siempre indica bounds con method='bounded' cuando la función no tiene un mínimo global (es decir, no está acotada):

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

Búsqueda de raíces

optimize.root_scalar encuentra dónde una función cruza el cero:

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

Interpolación

scipy.interpolate ajusta una curva suave a través de puntos de datos para que puedas estimar valores entre ellos.

Interpolación 1-D

interp1d crea una función de interpolación invocable a partir de pares discretos (x, y). El parámetro kind selecciona el método: 'linear' (predeterminado), '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)

La interpolación cúbica recupera el resultado exacto porque los datos provienen de un polinomio cuadrático, y el spline cúbico es lo suficientemente flexible para ajustarlo perfectamente.

Estadística

scipy.stats contiene más de 80 distribuciones de probabilidad continuas y discretas, además de una colección de pruebas estadísticas.

Estadística descriptiva y distribución normal

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

Pruebas de hipótesis

scipy.stats incluye pruebas t, pruebas chi-cuadrado, ANOVA, pruebas de Kolmogorov-Smirnov y más.

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

Variables aleatorias y muestreo

Cada distribución en scipy.stats expone la misma interfaz: pdf, cdf, ppf, rvs (variantes aleatorias) y 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}')

Procesamiento de Imágenes N-Dimensional

scipy.ndimage opera sobre arrays de cualquier dimensionalidad (imágenes, volúmenes, cubos de series temporales). A continuación se muestra un ejemplo autónomo que usa un array 2-D sintético para que no se necesite ningún archivo de imagen externo.

Desenfoque gaussiano y etiquetado de regiones conectadas

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

Operaciones comunes de `ndimage`

Función	Qué hace
`gaussian_filter(a, sigma)`	Suavizar con un kernel gaussiano
`sobel(a)`	Detectar bordes (gradiente de Sobel)
`label(a)`	Etiquetar regiones conectadas en un array binario
`binary_dilation(a)`	Expandir las regiones de primer plano
`zoom(a, factor)`	Redimensionar un array
`rotate(a, angle)`	Rotar un array (en grados)

Cuándo Usar SciPy vs. Otras Bibliotecas

Tarea	Herramienta preferida
Creación de arrays y matemáticas básicas	NumPy
DataFrames, series temporales, E/S	Pandas
Algoritmos científicos (integración, optimización)	SciPy
Aprendizaje automático	scikit-learn (se basa en SciPy)
Visualización	Matplotlib

SciPy no reemplaza a NumPy — depende de NumPy y lo extiende. En la práctica importarás ambos.

Referencia Rápida

Sub-paquete	Funciones principales
`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`