SciPy: Computação Científica Avançada em Python

O SciPy é a biblioteca de computação científica que estende o NumPy com algoritmos numéricos consagrados: otimização, ajuste de curvas, estatística, integração de equações diferenciais, interpolação, processamento de sinais e matrizes esparsas. Enquanto o NumPy fornece a estrutura de dados (o ndarray) e as operações elementares, o SciPy fornece os métodos numéricos de alto nível, em sua maioria apoiados em bibliotecas de referência em Fortran e C (LAPACK, MINPACK, FFTPACK). Este guia percorre os submódulos mais utilizados com exemplos completos e executáveis.

SciPy em comparação com o NumPy

A fronteira é simples de lembrar:

• NumPy: o array ndarray, a aritmética vetorizada, a álgebra linear básica, a geração de números aleatórios.
• SciPy: os algoritmos que operam sobre esses arrays, com implementações robustas fruto de décadas de computação científica.

O SciPy depende do NumPy e é importado submódulo por submódulo:

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

# Verificar a versão instalada
import scipy
print(scipy.__version__)     # 1.14.1 (por exemplo)

Existem algumas sobreposições (scipy.linalg estende numpy.linalg com decomposições adicionais; scipy.fft substitui com vantagem numpy.fft), e nesses casos a versão do SciPy costuma ser mais completa.

scipy.optimize: minimização e ajuste de curvas

Minimizar uma função com minimize

minimize encontra um mínimo local de uma função escalar de várias variáveis.

from scipy.optimize import minimize

# Função de Rosenbrock, um benchmark clássico
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] -> mínimo em (1, 1)
print(resultat.fun)          # 5.9e-10 (valor no mínimo, quase nulo)
print(resultat.nit)          # número de iterações

A escolha do método depende do problema: "Nelder-Mead" não requer gradiente, "BFGS" (padrão sem restrições) estima o gradiente por diferenças finitas, "L-BFGS-B" aceita limites e "SLSQP" lida com restrições gerais.

# Minimização com limites: x em [0, 2], y em [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 a raiz de uma equação em vez de um mínimo, use optimize.root_scalar ou optimize.brentq.

Ajustar um modelo a dados com curve_fit

curve_fit estima os parâmetros de um modelo por mínimos quadrados não lineares. Exemplo: ajustar um decaimento exponencial a medições ruidosas.

from scipy.optimize import curve_fit

# Modelo: a primeira variável é sempre x, as seguintes são os parâmetros
def modele(x, a, tau, c):
    return a * np.exp(-x / tau) + c

# Dados sintéticos ruidosos
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)

# Ajuste (p0: estimativa inicial dos parâmetros)
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

A diagonal da matriz de covariância fornece a variância de cada parâmetro estimado: sua raiz quadrada dá uma incerteza de um desvio padrão. Forneça sempre um p0 razoável para modelos não lineares, sob pena de convergir para um ótimo ruim.

scipy.stats: distribuições e testes estatísticos

Trabalhar com distribuições

Cada distribuição expõe a mesma interface: pdf (densidade), cdf (função de distribuição acumulada), ppf (quantis, inversa da cdf) e rvs (amostragem).

from scipy import stats

# Distribuição normal de média 100 e desvio padrão 15
loi = stats.norm(loc=100, scale=15)

print(loi.pdf(100))          # 0.02659... (densidade no pico)
print(loi.cdf(115))          # 0.8413... (P(X <= 115))
print(loi.ppf(0.975))        # 129.399... (quantil de 97,5 %)
print(loi.rvs(size=3, random_state=42))   # [107.45  84.44 109.72]

# Outras distribuições do mesmo molde: stats.t, stats.chi2, stats.binom,
# stats.poisson, stats.expon, stats.uniform...

Teste t: comparar duas médias

rng = np.random.default_rng(1)
groupe_a = rng.normal(loc=52, scale=8, size=40)   # ex.: notas com o método A
groupe_b = rng.normal(loc=48, scale=8, size=40)   # ex.: notas com o método 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 (por exemplo)

if p_valeur < 0.05:
    print("Différence significative au seuil de 5 %")
# Variante pareada (medições antes/depois nos mesmos sujeitos): stats.ttest_rel
# Variâncias desiguais: ttest_ind(..., equal_var=False) (teste de Welch)

Teste qui-quadrado: independência de duas variáveis categóricas

# Tabela de contingência: linhas = treinamento realizado, colunas = aprovação
observations = np.array([[45, 15],    # com treinamento: 45 aprovações, 15 reprovações
                         [30, 30]])   # sem treinamento: 30 aprovações, 30 reprovações

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)              # frequências esperadas sob a hipótese de independência
# [[37.5 22.5]
#  [37.5 22.5]]

Correlações

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

r, p = stats.pearsonr(x, y)          # correlação linear
print(f"Pearson  r = {r:.3f}, p = {p:.2e}")   # r = 0.812, p = 1.5e-24

rho, p = stats.spearmanr(x, y)       # correlação de postos (relações monotônicas)
print(f"Spearman rho = {rho:.3f}")   # rho = 0.795

Use Spearman quando a relação for monotônica, mas não linear, ou na presença de outliers.

scipy.integrate: integrais e equações diferenciais

Integração numérica com quad

from scipy.integrate import quad

# Integral de exp(-x²) entre 0 e o infinito (vale 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 (estimativa do erro)

Resolver uma EDO com solve_ivp

Exemplo: o oscilador harmônico amortecido, descrito pela equação de segunda ordem x” + 2ζω x’ + ω² x = 0, reescrita como um sistema de primeira ordem.

from scipy.integrate import solve_ivp

omega = 2.0    # frequência angular natural
zeta = 0.15    # razão de amortecimento

def oscillateur(t, y):
    x, v = y                              # posição, velocidade
    return [v, -2 * zeta * omega * v - omega ** 2 * x]

solution = solve_ivp(
    oscillateur,
    t_span=(0, 20),                       # intervalo de integração
    y0=[1.0, 0.0],                        # posição 1, velocidade nula em t=0
    t_eval=np.linspace(0, 20, 400),       # pontos onde avaliar a solução
    method="RK45",                        # Runge-Kutta adaptativo (padrão)
)

print(solution.success)                   # True
print(solution.y.shape)                   # (2, 400) -> posição e velocidade
print(solution.y[0, :3])                  # [1.  0.995 0.980] (posição amortecida)

# Gráfico da posição ao longo do tempo
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 sistemas rígidos (constantes de tempo muito diferentes), use method="Radau" ou method="BDF".

scipy.interpolate: interpolação de dados

A interpolação constrói uma função contínua que passa por pontos de medição discretos.

from scipy.interpolate import CubicSpline, interp1d

# Pontos de medição esparsos
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 cúbica: suave, derivável, recomendada
spline = CubicSpline(x_mesures, y_mesures)

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

print(spline(2.5))                   # 0.5546... (valor interpolado entre 2 e 3)
print(spline(2.5, 1))                # derivada primeira no mesmo ponto

# Interpolação linear simples
lineaire = interp1d(x_mesures, y_mesures, kind="linear")
print(lineaire(2.5))                 # 0.5 (média dos pontos vizinhos)

Atenção: a interpolação só é confiável dentro do domínio das medições. Para extrapolar ou suavizar dados ruidosos, prefira um ajuste de modelo (curve_fit) ou uma spline de suavização (scipy.interpolate.make_smoothing_spline).

scipy.signal e scipy.fft: processamento de sinais

Análise espectral com a FFT

from scipy.fft import fft, fftfreq

# Sinal amostrado: soma de duas senoides (5 Hz e 50 Hz) mais ruído
fe = 500                              # frequência de amostragem (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))

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

# Amplitudes nas frequências positivas
moitie = t.size // 2
amplitudes = 2 / t.size * np.abs(spectre[:moitie])

pics = frequences[:moitie][amplitudes > 0.3]
print(pics)                           # [ 5. 50.] -> as duas frequências recuperadas

Filtragem com um filtro de Butterworth

from scipy.signal import butter, filtfilt

# Filtro passa-baixas de ordem 4, frequência de corte de 20 Hz
b, a = butter(N=4, Wn=20, btype="low", fs=fe)

# filtfilt aplica o filtro nos dois sentidos: nenhuma defasagem
signal_filtre = filtfilt(b, a, signal_brut)

# A componente de 50 Hz e o ruído de alta frequência são eliminados,
# a senoide de 5 Hz é preservada.

O scipy.signal também oferece detecção de picos (find_peaks), espectrogramas (ShortTimeFFT), convolução e o projeto de diversos 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 -> um pico por período da senoide de 5 Hz

scipy.sparse: matrizes esparsas

Quando uma matriz contém uma esmagadora maioria de zeros (grafos, sistemas de equações oriundos de malhas, matrizes termo-documento), armazená-la de forma densa desperdiça memória e processamento.

from scipy import sparse

# Matriz esparsa no formato COO (construção) e depois CSR (cálculo)
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 elementos não nulos armazenados (de 9)

# Produto matriz-vetor: apenas os elementos não nulos são percorridos
v = np.array([1.0, 2.0, 3.0])
print(m @ v)                          # [ 6. 12. 23.]

# Resolução de um sistema linear esparso
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.]

Ordens de grandeza: uma matriz densa de 100.000 × 100.000 em float64 ocuparia 80 GB; a mesma matriz com 0,01 % de elementos não nulos cabe em alguns megabytes no formato CSR.

Tabela resumo dos submódulos

Submódulo	Domínio	Funções principais
scipy.optimize	Otimização, ajuste, raízes	minimize, curve_fit, root_scalar, linprog
scipy.stats	Distribuições, testes, correlações	norm, ttest_ind, chi2_contingency, pearsonr
scipy.integrate	Integrais, equações diferenciais	quad, solve_ivp
scipy.interpolate	Interpolação, splines	CubicSpline, interp1d, make_smoothing_spline
scipy.signal	Filtragem, picos, espectrogramas	butter, filtfilt, find_peaks
scipy.fft	Transformadas de Fourier	fft, ifft, fftfreq, rfft
scipy.sparse	Matrizes esparsas	csr_array, coo_array, sparse.linalg.spsolve
scipy.linalg	Álgebra linear estendida	lu, svd, expm, solve_banded
scipy.spatial	Geometria, distâncias	KDTree, distance.cdist, ConvexHull
scipy.ndimage	Processamento de imagens n-D	gaussian_filter, label, rotate

Boas práticas

• Importe os submódulos explicitamente (from scipy import optimize): apenas import scipy não carrega os submódulos.
• Sempre verifique resultat.success após minimize ou solve_ivp: um algoritmo que não convergiu retorna um resultado mesmo assim.
• Forneça uma estimativa inicial p0 sensata ao curve_fit: os mínimos quadrados não lineares convergem para um ótimo local que depende do ponto de partida.
• Interprete os valores-p com cautela: um valor-p abaixo de 0,05 não mede nem o tamanho do efeito nem sua importância prática; relate também as médias e os intervalos de confiança.
• Escolha o método adequado ao problema: Radau ou BDF para EDOs rígidas, L-BFGS-B sempre que houver limites, o teste de Welch quando as variâncias diferem.
• Não extrapole com uma interpolação: fora do domínio das medições, uma spline cúbica diverge rapidamente.
• Migre para matrizes esparsas assim que a densidade de elementos não nulos cair abaixo de alguns por cento em matrizes grandes.

Conclusão

O SciPy transforma o Python em um ambiente completo de computação científica: otimização e ajuste de modelos com optimize, inferência estatística com stats, resolução de equações diferenciais com integrate, reconstrução de funções com interpolate, análise de frequência e filtragem com signal e fft, e escalabilidade com sparse. Cada submódulo encapsula algoritmos de referência validados há décadas, o que evita que você reimplemente métodos numéricos delicados. Combinado ao NumPy para as estruturas de dados, ao Pandas para tabelas rotuladas e ao Matplotlib para visualização, o SciPy completa o conjunto de ferramentas padrão de qualquer trabalho científico ou de engenharia em Python.