SciPy-Tutorial — Wissenschaftliches Rechnen in Python

SciPy (Scientific Python) ist eine Open-Source-Bibliothek, die auf NumPy aufbaut und eine umfangreiche Sammlung von Algorithmen für Mathematik, Wissenschaft und Ingenieurwesen bereitstellt. Während NumPy das ndarray und grundlegende Operationen liefert, bietet SciPy spezialisierte Löser: numerische Integration, Optimierung, Interpolation, Signalverarbeitung, Statistik, räumliche Algorithmen und mehr.

Dieses Kapitel behandelt:

Installation von SciPy und die Import-Konvention
Lineare Algebra (scipy.linalg)
Numerische Integration und Differentiation (scipy.integrate)
Optimierung — Minima und Nullstellen finden (scipy.optimize)
Interpolation (scipy.interpolate)
Statistik (scipy.stats)
N-dimensionale Bildverarbeitung (scipy.ndimage)

SciPy installieren

SciPy ist in der Anaconda-Distribution enthalten. So installieren Sie es mit pip:

pip install scipy

SciPy hängt von NumPy ab, das pip automatisch installiert, sofern es noch nicht vorhanden ist.

Import-Konvention

SciPy ist in Teilpakete unterteilt. Importieren Sie nur die benötigten Teilpakete anstatt die gesamte Bibliothek:

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

Sie können auch die installierte Version abfragen:

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

Lineare Algebra

scipy.linalg erweitert die linearen Algebra-Routinen von NumPy um zusätzliche Zerlegungen und Löser. Alle Funktionen arbeiten mit regulären NumPy-Arrays.

Determinante und Inverse

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. Die Inverse erfüllt a @ inv == I.

Eigenwerte und Eigenvektoren

Eigenwerte beschreiben, wie eine Matrix den Raum streckt; Eigenvektoren geben die Richtungen an, die sich nicht drehen.

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

Jede Spalte von eigenvectors entspricht dem zugehörigen Eigenwert. Die Eigenwerte werden als komplexe Zahlen zurückgegeben, auch wenn der Imaginärteil null ist.

Singulärwertzerlegung (SVD)

Die SVD zerlegt eine Matrix A in drei Matrizen U, s, Vt, sodass A = U @ diag(s) @ Vt gilt. Sie ist die Grundlage der Hauptkomponentenanalyse (PCA) und vieler Dimensionsreduktionsverfahren.

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

Ein lineares Gleichungssystem lösen

linalg.solve ist die korrekte Methode, um Ax = b zu lösen. Sie ist schneller und numerisch stabiler als die Berechnung der Inversen mit anschließender Multiplikation.

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

Numerische Integration

scipy.integrate stellt Routinen zur Berechnung bestimmter Integrale bereit, wenn eine analytische Lösung unpraktisch ist.

Einvariablen-Integration mit `quad`

integrate.quad verwendet adaptive Quadratur, um eine Funktion über ein Intervall zu integrieren. Es gibt das Ergebnis und eine Schätzung des absoluten Fehlers zurück.

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)

Numerische Differentiation

approx_fprime von SciPy berechnet einen Finite-Differenzen-Gradienten. Für skalare Funktionen ist die zentrale Differenz derivative aus scipy.misc einfacher:

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

Optimierung

scipy.optimize findet Minima, Maxima (durch Minimierung des Negativen) und Nullstellen von Funktionen.

Minimierung einer multivariaten Funktion

optimize.minimize unterstützt viele Methoden (Nelder-Mead, BFGS, L-BFGS-B, …). Die Standardmethode wird automatisch gewählt.

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)

Minimierung einer skalaren Funktion über ein begrenztes Intervall

optimize.minimize_scalar ist einfacher für Einvariablenprobleme. Geben Sie immer bounds mit method='bounded' an, wenn die Funktion kein globales Minimum hat (d. h. unbeschränkt ist):

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

Nullstellen finden

optimize.root_scalar findet, wo eine Funktion die Nulllinie schneidet:

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

Interpolation

scipy.interpolate legt eine glatte Kurve durch Datenpunkte, sodass Sie Werte zwischen ihnen schätzen können.

1-D-Interpolation

interp1d erzeugt eine aufrufbare Interpolationsfunktion aus diskreten (x, y)-Paaren. Der Parameter kind wählt die Methode: 'linear' (Standard), 'quadratic' oder '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)

Die kubische Interpolation liefert das exakte Ergebnis, weil die Daten aus einem quadratischen Polynom stammen und der kubische Spline flexibel genug ist, dieses perfekt anzupassen.

Statistik

scipy.stats enthält über 80 kontinuierliche und diskrete Wahrscheinlichkeitsverteilungen sowie eine Sammlung statistischer Tests.

Deskriptive Statistik und Normalverteilung

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

Hypothesentests

scipy.stats enthält t-Tests, Chi-Quadrat-Tests, ANOVA, Kolmogorov-Smirnov-Tests und mehr.

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

Zufallsvariablen und Stichproben

Jede Verteilung in scipy.stats bietet dieselbe Schnittstelle: pdf, cdf, ppf, rvs (Zufallsstichproben) und 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}')

N-dimensionale Bildverarbeitung

scipy.ndimage arbeitet mit Arrays beliebiger Dimensionalität (Bilder, Volumina, Zeit-Reihen-Würfel). Das folgende in sich abgeschlossene Beispiel verwendet ein synthetisches 2-D-Array, sodass keine externe Bilddatei benötigt wird.

Gaußsche Weichzeichnung und Beschriftung verbundener Regionen

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

Häufige `ndimage`-Operationen

Funktion	Beschreibung
`gaussian_filter(a, sigma)`	Weichzeichnung mit einem Gaußschen Kern
`sobel(a)`	Kantenerkennung (Sobel-Gradient)
`label(a)`	Beschriftung verbundener Regionen in einem binären Array
`binary_dilation(a)`	Vordergrundregionen ausdehnen
`zoom(a, factor)`	Array skalieren
`rotate(a, angle)`	Array drehen (in Grad)

Wann SciPy und wann andere Bibliotheken verwenden

Aufgabe	Bevorzugtes Werkzeug
Array-Erstellung und grundlegende Mathematik	NumPy
DataFrames, Zeitreihen, Ein-/Ausgabe	Pandas
Wissenschaftliche Algorithmen (Integration, Optimierung)	SciPy
Maschinelles Lernen	scikit-learn (basiert auf SciPy)
Visualisierung	Matplotlib

SciPy ist kein Ersatz für NumPy — es hängt von NumPy ab und erweitert es. In der Praxis werden Sie beide importieren.

Kurzreferenz

Teilpaket	Wichtige Funktionen
`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`