この記事の対象読者
- NumPyは使えるが、SciPyとの違いがわからない方
- 最適化・統計検定・信号処理・補間など「数学寄りの処理」が必要になった方
-
scipy.optimizeやscipy.statsを使ったことはあるが、全体像が掴めていない方 - 科学技術計算や研究開発でPythonを使い始めた方
この記事で得られること
この記事を読むと、以下のことが理解できます:
- SciPyとNumPyの「棲み分け」— なぜ両方必要なのか
- 主要サブモジュールの「役割と使いどころ」— optimize / stats / signal / interpolate / linalg / integrate / sparse
- 最適化の「仕組み」— minimize関数が内部でやっていること
- 統計検定の「使い分け」— t検定 / カイ二乗検定 / 正規性検定をいつ使うか
- 疎行列の「威力」— 大規模データでメモリを99%削減するテクニック
この記事で扱わないこと
- 各関数の数学的導出の詳細
- SymPy による記号計算
- MATLAB / Mathematica との網羅的な比較
本記事ではSciPy 1.12〜1.14系を前提としています。コード例はPython 3.11以降で動作確認済みです。
本記事の比喩について
この記事では、SciPyの仕組みを「大学の理工学部棟」に見立てて解説する。
| SciPyの概念 | 理工学部棟での対応物 |
|---|---|
| SciPy全体 | 理工学部棟 — 各フロアに専門の研究室が入っている |
| NumPy | 学部棟の共通インフラ — 電気・水道・ネットワーク。全研究室が依存 |
scipy.optimize |
最適化研究室 — 最もコスパの良い解を見つける専門家 |
scipy.stats |
統計学研究室 — データの分布を分析し、仮説を検証する |
scipy.signal |
信号処理研究室 — 波形データを解析・フィルタリングする |
scipy.interpolate |
補間研究室 — データの隙間を滑らかに埋める |
scipy.linalg |
線形代数研究室 — 行列演算の専門家集団 |
scipy.integrate |
積分研究室 — 面積や体積を計算する |
scipy.sparse |
疎行列研究室 — ほとんどゼロのデータを効率的に扱う |
各研究室はそれぞれ独立した専門分野を持つが、共通インフラとしてNumPyの配列を使っている。この構造を意識しながら読み進めてほしい。
1. SciPyとNumPyの「棲み分け」— なぜ両方必要なのか
1.1 NumPyは「インフラ」、SciPyは「専門研究室」
NumPyは多次元配列と基本的な数学演算を提供する。四則演算、行列積、線形代数の基本、乱数生成。これだけで多くの計算は可能だ。
ではSciPyは何をするのか。理工学部棟の比喩に戻ろう。NumPyが提供するのは電気・水道・ネットワークといった共通インフラだ。どの研究室もこのインフラの上で動いている。SciPyは、その上に構築された専門の研究設備を提供する。
1.2 具体例で見る棲み分け
| やりたいこと | NumPyでできる？ | SciPyが必要？ |
|---|---|---|
| 配列の四則演算・行列積 | ✅ | 不要 |
| 平均・分散・標準偏差 | ✅ | 不要 |
| 固有値分解・逆行列 | ✅ 基本的なもの | ✅ 高度なもの |
| 非線形方程式の求解 | ❌ | ✅ optimize |
| 統計検定・分布フィッティング | ❌ | ✅ stats |
| 信号のフィルタリング・FFT | ❌ | ✅ signal |
| データの補間・スプライン | ❌ | ✅ interpolate |
| 数値積分 | ❌ | ✅ integrate |
| 疎行列の効率的な演算 | ❌ | ✅ sparse |
NumPyにも numpy.linalg や numpy.fft がありますが、SciPyの対応モジュールはより高機能・高速な実装を含んでいます。迷ったらSciPy側を使うのが安全です。
SciPyの立ち位置が明確になったところで、各研究室を順番に訪問していこう。最初は実務で最も頻出する「最適化研究室」からだ。
2. scipy.optimize — 最適化研究室
2.1 最適化とは何か
最適化とは、ある関数の値を最小または最大にする入力を見つける問題だ。
理工学部棟の最適化研究室は、あらゆる「ベストな答え探し」を引き受ける。製造コストを最小化する生産計画、投資リターンを最大化するポートフォリオ、実験データに最もフィットする曲線のパラメータ。すべて最適化問題だ。
2.2 minimize — 関数の最小値を見つける
from scipy.optimize import minimize
import numpy as np
# 最小化したい関数（ローゼンブロック関数 — 最適化のベンチマーク関数）
def rosenbrock(x):
return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2
# 初期値から最適化を開始
x0 = np.array([-1.0, 1.0])
result = minimize(rosenbrock, x0, method="Nelder-Mead")
print(f"最適解: {result.x}") # [1.0, 1.0] に近い値
print(f"最小値: {result.fun}") # 0 に近い値
print(f"成功: {result.success}") # True
print(f"反復回数: {result.nit}")
2.3 最適化手法の使い分け
| 手法 | method= |
特徴 | 使いどころ |
|---|---|---|---|
| Nelder-Mead | "Nelder-Mead" |
勾配不要、ロバスト | 微分できない関数、まず試す手法 |
| BFGS | "BFGS" |
準ニュートン法、高速 | 滑らかな関数の高速最適化 |
| L-BFGS-B | "L-BFGS-B" |
メモリ効率が良い、境界条件対応 | 高次元問題、変数に範囲制約あり |
| Powell | "Powell" |
勾配不要 | 非滑らかな関数 |
| trust-constr | "trust-constr" |
等式/不等式制約対応 | 制約付き最適化問題 |
2.4 curve_fit — 実験データへの曲線フィッティング
科学計算で最も使われるSciPyの機能の一つが curve_fit だ。
from scipy.optimize import curve_fit
import numpy as np
# 実験データ（ノイズ付き指数減衰）
np.random.seed(42)
t = np.linspace(0, 5, 50)
y_true = 3.0 * np.exp(-1.5 * t) + 0.5
y_data = y_true + 0.1 * np.random.randn(len(t))
# フィットしたいモデル関数
def exp_decay(t, a, b, c):
"""y = a * exp(-b * t) + c"""
return a * np.exp(-b * t) + c
# カーブフィッティング
popt, pcov = curve_fit(exp_decay, t, y_data)
a_fit, b_fit, c_fit = popt
perr = np.sqrt(np.diag(pcov)) # パラメータの標準誤差
print(f"a = {a_fit:.3f} ± {perr[0]:.3f} (真値: 3.0)")
print(f"b = {b_fit:.3f} ± {perr[1]:.3f} (真値: 1.5)")
print(f"c = {c_fit:.3f} ± {perr[2]:.3f} (真値: 0.5)")
curve_fit は初期値に敏感です。フィッティングが収束しない場合は p0 パラメータで妥当な初期値を与えてください。物理的な意味から初期値を推定するのが定石です。
最適化研究室を出て、次は統計学研究室を訪問しよう。「このデータの差は偶然か、それとも意味があるのか」を判定する統計検定の世界だ。
3. scipy.stats — 統計学研究室
3.1 確率分布 — データの「形」を知る
scipy.stats には100以上の確率分布が実装されている。データがどんな「形」をしているかを知ることは、適切な分析手法を選ぶ第一歩だ。
from scipy import stats
import numpy as np
# 正規分布 N(μ=170, σ=6) — 日本人成人男性の身長分布に近似
height_dist = stats.norm(loc=170, scale=6)
# 確率密度関数の値
print(f"身長170cmの確率密度: {height_dist.pdf(170):.4f}")
# 累積分布関数 — 「180cm以下の割合は？」
print(f"180cm以下の割合: {height_dist.cdf(180):.4f}") # ≈ 0.952
# パーセンタイル — 「上位5%は何cm以上？」
print(f"上位5%の境界: {height_dist.ppf(0.95):.1f}cm") # ≈ 179.9cm
# ランダムサンプリング
samples = height_dist.rvs(size=1000, random_state=42)
print(f"サンプル平均: {samples.mean():.1f}cm, 標準偏差: {samples.std():.1f}cm")
3.2 統計検定の使い分け
| 検定 | 関数 | 問い | 使いどころ |
|---|---|---|---|
| t検定 | stats.ttest_ind |
2群の平均に差があるか？ | A/Bテスト、薬効比較 |
| 対応のあるt検定 | stats.ttest_rel |
同じ対象の前後比較で差があるか？ | 施策前後の効果測定 |
| カイ二乗検定 | stats.chi2_contingency |
2つのカテゴリ変数に関連があるか？ | アンケート分析 |
| 正規性検定 | stats.shapiro |
データは正規分布に従うか？ | t検定の前提確認 |
| マン・ホイットニーU検定 | stats.mannwhitneyu |
正規分布でない2群の差は？ | 正規性が仮定できない場合 |
| KS検定 | stats.ks_2samp |
2つの分布は同じか？ | 分布の一致判定 |
3.3 実践: A/Bテストの統計検定
from scipy import stats
import numpy as np
np.random.seed(42)
# Webサイトの2つのデザインのコンバージョン率
# デザインA: 1000人中120人がコンバージョン
# デザインB: 1000人中145人がコンバージョン
a_conversions = np.array([1]*120 + [0]*880)
b_conversions = np.array([1]*145 + [0]*855)
# 正規性の確認（大標本なのでCLTで正規近似可能だが念のため）
# t検定の実行
t_stat, p_value = stats.ttest_ind(a_conversions, b_conversions)
print(f"t統計量: {t_stat:.4f}")
print(f"p値: {p_value:.4f}")
if p_value < 0.05:
print("→ 有意水準5%で有意差あり。デザインBの方がコンバージョン率が高い。")
else:
print("→ 有意差なし。差は偶然の範囲。")
# 効果量（Cohen's d）も確認
mean_diff = b_conversions.mean() - a_conversions.mean()
pooled_std = np.sqrt((a_conversions.std()**2 + b_conversions.std()**2) / 2)
cohens_d = mean_diff / pooled_std
print(f"効果量 (Cohen's d): {cohens_d:.4f}")
p値だけで判断するのは危険です。効果量と信頼区間も併せて確認してください。「統計的に有意だが効果量が極めて小さい」ケースは実務では意味がないことが多いです。
3.4 分布フィッティング — データの分布を推定する
from scipy import stats
import numpy as np
# 未知の分布に従うデータ（実は対数正規分布）
np.random.seed(42)
data = stats.lognorm.rvs(s=0.5, loc=0, scale=np.exp(3), size=1000)
# 候補分布にフィッティングして最も当てはまりの良いものを探す
distributions = [stats.norm, stats.lognorm, stats.expon, stats.gamma]
results = []
for dist in distributions:
try:
params = dist.fit(data)
# KS検定でフィットの良さを評価
ks_stat, ks_p = stats.kstest(data, dist.cdf, args=params)
results.append({
"distribution": dist.name,
"ks_statistic": ks_stat,
"ks_p_value": ks_p,
"params": params
})
except Exception:
pass
# p値が最も大きい分布が最もフィットが良い
results.sort(key=lambda x: x["ks_p_value"], reverse=True)
best = results[0]
print(f"最もフィットの良い分布: {best['distribution']}")
print(f"KS統計量: {best['ks_statistic']:.4f}, p値: {best['ks_p_value']:.4f}")
統計学研究室を出て、次は信号処理研究室を訪問しよう。時系列データやセンサーデータを扱う際に必須となるフィルタリングとFFTの世界だ。
4. scipy.signal — 信号処理研究室
4.1 フィルタリング — ノイズを除去する
センサーデータや音声データには必ずノイズが混じる。信号処理研究室の仕事は、このノイズを除去して本来の信号を取り出すことだ。
from scipy import signal
import numpy as np
# サンプル信号: 5Hzの正弦波 + ノイズ
np.random.seed(42)
fs = 1000 # サンプリング周波数 (Hz)
t = np.arange(0, 1, 1/fs)
clean_signal = np.sin(2 * np.pi * 5 * t) # 5Hzの信号
noisy_signal = clean_signal + 0.5 * np.random.randn(len(t)) # ノイズ追加
# バターワースローパスフィルタ（10Hz以上をカット）
b, a = signal.butter(N=4, Wn=10, fs=fs, btype="low")
filtered = signal.filtfilt(b, a, noisy_signal)
print(f"ノイズ付き信号のRMS: {np.sqrt(np.mean(noisy_signal**2)):.4f}")
print(f"フィルタ後のRMS: {np.sqrt(np.mean(filtered**2)):.4f}")
print(f"元信号のRMS: {np.sqrt(np.mean(clean_signal**2)):.4f}")
4.2 FFT — 周波数成分を分析する
from scipy import fft
import numpy as np
# 複数の周波数成分を持つ信号
fs = 1000
t = np.arange(0, 1, 1/fs)
signal_data = (
1.0 * np.sin(2 * np.pi * 50 * t) + # 50Hz成分（振幅1.0）
0.5 * np.sin(2 * np.pi * 120 * t) + # 120Hz成分（振幅0.5）
0.3 * np.sin(2 * np.pi * 300 * t) # 300Hz成分（振幅0.3）
)
# FFTの実行
N = len(signal_data)
yf = fft.fft(signal_data)
xf = fft.fftfreq(N, 1/fs)
# 正の周波数のみ取得し、振幅スペクトルを計算
mask = xf >= 0
freqs = xf[mask]
amplitudes = 2.0 / N * np.abs(yf[mask])
# ピーク検出
peaks, _ = signal.find_peaks(amplitudes, height=0.1)
for p in peaks:
print(f" 周波数: {freqs[p]:.0f}Hz, 振幅: {amplitudes[p]:.2f}")
scipy.fft はNumPyの numpy.fft より高速な実装です。大規模データのFFTでは scipy.fft を使ってください。
5. scipy.interpolate — 補間研究室
5.1 補間とは何か
補間は既知のデータ点の間を滑らかに埋める操作だ。理工学部棟の補間研究室は、測定点が離散的なデータから連続的な曲線を復元する専門家だ。
from scipy.interpolate import interp1d, CubicSpline
import numpy as np
# 離散的な測定データ（6時間おきの気温）
hours = np.array([0, 6, 12, 18, 24])
temps = np.array([15, 18, 28, 22, 16])
# 線形補間
linear_interp = interp1d(hours, temps, kind="linear")
# 3次スプライン補間（滑らか）
cubic_interp = CubicSpline(hours, temps)
# 1時間刻みで補間値を取得
hours_fine = np.arange(0, 24.1, 1)
temps_linear = linear_interp(hours_fine)
temps_cubic = cubic_interp(hours_fine)
print("時刻 線形補間 スプライン補間")
for h in [3, 9, 15, 21]:
print(f" {h:2d}時 {linear_interp(h):.1f}℃ {cubic_interp(h):.1f}℃")
5.2 補間手法の使い分け
| 手法 | 関数 | 特徴 | 使いどころ |
|---|---|---|---|
| 線形補間 | interp1d(kind="linear") |
直線で結ぶ。単純で高速 | ざっくりした補間 |
| 3次スプライン | CubicSpline |
滑らかな曲線。微分も可能 | 物理データの補間 |
| Akima補間 | Akima1DInterpolator |
振動が少ない | 振動を避けたいとき |
| 2次元補間 | RegularGridInterpolator |
2D格子データの補間 | 地理データ、画像のリサンプリング |
6. scipy.linalg — 線形代数研究室
6.1 NumPyのlinalgとの違い
numpy.linalg でも基本的な線形代数はできるが、scipy.linalg はより高度で高速な実装を含んでいる。
from scipy import linalg
import numpy as np
A = np.array([[1, 2], [3, 4]])
b = np.array([5, 6])
# 連立方程式 Ax = b の求解
x = linalg.solve(A, b)
print(f"解: {x}") # [-4. 4.5]
# 検算: A @ x = b
print(f"検算: {A @ x}") # [5. 6.]
# LU分解
P, L, U = linalg.lu(A)
print(f"L:\n{L}")
print(f"U:\n{U}")
# 特異値分解 (SVD)
U_svd, s, Vt = linalg.svd(A)
print(f"特異値: {s}")
6.2 scipy.linalgが必要になる場面
| 操作 | numpy.linalg | scipy.linalg |
|---|---|---|
| 基本の固有値分解 | ✅ eig
|
✅ eig
|
| 帯行列・三角行列の特殊解法 | ❌ | ✅ solve_banded, solve_triangular
|
| コレスキー分解 | ✅ cholesky
|
✅ cholesky + α |
| 行列指数関数・対数関数 | ❌ | ✅ expm, logm
|
| 行列方程式 | ❌ | ✅ solve_sylvester, solve_lyapunov
|
「迷ったら scipy.linalg を使う」が安全な方針です。NumPy版よりも常に同等以上の機能と性能を提供します。
7. scipy.integrate — 積分研究室
7.1 数値積分
解析的に積分できない関数でも、数値的に積分値を計算できる。
from scipy import integrate
import numpy as np
# 例: ガウス関数の積分 ∫₋∞^∞ e^(-x²) dx = √π
def gaussian(x):
return np.exp(-x**2)
result, error = integrate.quad(gaussian, -np.inf, np.inf)
print(f"積分値: {result:.6f}")
print(f"理論値: {np.sqrt(np.pi):.6f}")
print(f"誤差: {error:.2e}")
7.2 常微分方程式の数値解法
from scipy.integrate import solve_ivp
import numpy as np
# 減衰振動: d²x/dt² + 2γ(dx/dt) + ω²x = 0
# → 1階連立: dx/dt = v, dv/dt = -2γv - ω²x
def damped_oscillator(t, state, gamma, omega):
x, v = state
dxdt = v
dvdt = -2 * gamma * v - omega**2 * x
return [dxdt, dvdt]
gamma = 0.1 # 減衰係数
omega = 2.0 # 角周波数
sol = solve_ivp(
damped_oscillator,
t_span=[0, 20],
y0=[1.0, 0.0], # 初期位置1、初速0
args=(gamma, omega),
t_eval=np.linspace(0, 20, 500),
method="RK45"
)
print(f"計算点数: {len(sol.t)}")
print(f"最終位置: {sol.y[0, -1]:.6f}")
print(f"最終速度: {sol.y[1, -1]:.6f}")
8. scipy.sparse — 疎行列研究室
8.1 疎行列とは何か
大規模な行列の中には、ほとんどの要素がゼロであるものが多い。SNSのフォロー関係、Webのリンク構造、推薦システムのユーザー×アイテム行列。これらはすべて疎行列だ。
理工学部棟の疎行列研究室は、「ゼロを保存しない」というシンプルなアイデアでメモリと計算時間を劇的に削減する。
from scipy import sparse
import numpy as np
# 通常の密行列: 10000×10000 = 1億要素
# → float64で約800MB
dense_size_mb = 10000 * 10000 * 8 / (1024**2)
print(f"密行列のメモリ: {dense_size_mb:.0f} MB")
# 疎行列: 非ゼロ要素が0.1%しかない場合
n_nonzero = 10000 * 10000 * 0.001 # 10万要素
sparse_size_mb = n_nonzero * (8 + 4 + 4) / (1024**2) # 値+行+列のインデックス
print(f"疎行列のメモリ: {sparse_size_mb:.1f} MB")
print(f"削減率: {(1 - sparse_size_mb / dense_size_mb) * 100:.1f}%")
8.2 疎行列の形式
| 形式 | クラス | 特徴 | 使いどころ |
|---|---|---|---|
| CSR | csr_matrix |
行方向アクセスが高速 | 行列×ベクトル積、行スライス |
| CSC | csc_matrix |
列方向アクセスが高速 | 列スライス |
| COO | coo_matrix |
構築が高速 | 行列の初期構築 |
| LIL | lil_matrix |
要素の追加が高速 | 逐次的な行列構築 |
from scipy import sparse
import numpy as np
# COO形式で疎行列を構築（非ゼロ要素だけ指定）
row = [0, 0, 1, 2, 2, 3]
col = [0, 2, 1, 0, 2, 3]
data = [1, 2, 3, 4, 5, 6]
coo = sparse.coo_matrix((data, (row, col)), shape=(4, 4))
print(f"形式: {type(coo).__name__}")
print(f"非ゼロ要素数: {coo.nnz}")
print(f"密度: {coo.nnz / (4*4):.2%}")
# 演算にはCSR形式に変換
csr = coo.tocsr()
result = csr @ np.array([1, 1, 1, 1]) # 行列×ベクトル積
print(f"行列×ベクトル積: {result}")
# 密行列との相互変換
dense = csr.toarray()
print(f"密行列:\n{dense}")
8.3 疎行列の実用例: 推薦システムの類似度計算
from scipy import sparse
from sklearn.metrics.pairwise import cosine_similarity
import numpy as np
# ユーザー×映画の評価行列（10000ユーザー × 5000映画、密度1%）
np.random.seed(42)
n_users, n_movies = 10000, 5000
density = 0.01
# 疎行列で生成
ratings = sparse.random(n_users, n_movies, density=density, format="csr", random_state=42)
ratings.data = np.random.randint(1, 6, size=ratings.nnz).astype(float) # 1-5の評価
print(f"行列サイズ: {ratings.shape}")
print(f"非ゼロ要素: {ratings.nnz:,}")
print(f"密度: {ratings.nnz / (n_users * n_movies):.2%}")
print(f"疎行列メモリ: {(ratings.data.nbytes + ratings.indices.nbytes + ratings.indptr.nbytes) / 1024**2:.1f} MB")
print(f"密行列なら: {n_users * n_movies * 8 / 1024**2:.0f} MB")
疎行列を .toarray() で密行列に変換するとメモリが爆発するケースがあります。大規模データでは疎行列のまま演算することを徹底してください。
9. よくあるエラーと対処法
| # | エラー/症状 | 原因 | 対処法 |
|---|---|---|---|
| 1 |
optimize.minimize が収束しない |
初期値が悪い、または関数が非凸 | 初期値を変える、method を変える、tol を緩める |
| 2 |
curve_fit で OptimizeWarning: Covariance could not be estimated
|
パラメータが多すぎる、またはデータ点が少なすぎる |
p0 で適切な初期値を指定、パラメータ数を減らす |
| 3 |
stats.ttest_ind のp値が極端に小さい |
サンプルサイズが大きすぎて微小な差でも有意になる | 効果量を併せて確認。実質的な意味があるか判断する |
| 4 |
sparse.csr_matrix の要素追加が遅い |
CSR形式は構造変更に弱い | LIL形式で構築してからCSRに変換 |
| 5 |
linalg.solve で LinAlgError: singular matrix
|
行列が特異（逆行列が存在しない） |
linalg.lstsq で最小二乗解を求める |
| 6 | FFTの結果が期待と異なる | サンプリング周波数の指定漏れ |
fft.fftfreq(N, 1/fs) で正しい周波数軸を生成 |
| 7 |
integrate.quad の精度が悪い |
被積分関数が特異点を持つ |
points パラメータで特異点を指定、分割統治する |
10. SciPy環境診断スクリプト
#!/usr/bin/env python3
"""
SciPy環境診断スクリプト v1.0
- バージョン確認
- 各サブモジュールの動作テスト
- BLAS/LAPACKバックエンド確認
"""
import sys
import time
def diagnose_scipy():
try:
import scipy
import numpy as np
except ImportError as e:
print(f"❌ {e}")
print(" pip install scipy numpy でインストールしてください")
return
print("=" * 60)
print(" SciPy 環境診断レポート")
print("=" * 60)
# 1. バージョン
print("\n[1] バージョン情報")
print(f" SciPy: {scipy.__version__}")
print(f" NumPy: {np.__version__}")
print(f" Python: {sys.version.split()[0]}")
# 2. BLAS/LAPACKバックエンド
print("\n[2] 線形代数バックエンド")
try:
config = np.show_config(mode="dicts")
if isinstance(config, dict):
blas = config.get("Build Dependencies", {}).get("blas", {})
print(f" BLAS: {blas.get('name', '不明')} {blas.get('version', '')}")
else:
print(" （np.show_config()で確認）")
except Exception:
print(" （確認不可）")
# 3. 各サブモジュール動作テスト
print("\n[3] サブモジュール動作テスト")
tests = {
"optimize": _test_optimize,
"stats": _test_stats,
"signal": _test_signal,
"interpolate": _test_interpolate,
"linalg": _test_linalg,
"integrate": _test_integrate,
"sparse": _test_sparse,
}
for name, test_fn in tests.items():
try:
start = time.perf_counter()
result = test_fn()
elapsed = time.perf_counter() - start
print(f" ✅ scipy.{name:15s} {elapsed:.4f}秒 {result}")
except Exception as e:
print(f" ❌ scipy.{name:15s} エラー: {e}")
print("\n" + "=" * 60)
print(" 診断完了")
print("=" * 60)
def _test_optimize():
from scipy.optimize import minimize
import numpy as np
result = minimize(lambda x: (x[0]-1)**2 + (x[1]-2)**2, [0, 0])
return f"最適解: [{result.x[0]:.1f}, {result.x[1]:.1f}]"
def _test_stats():
from scipy import stats
import numpy as np
data = np.random.randn(1000)
stat, p = stats.normaltest(data)
return f"正規性検定 p={p:.4f}"
def _test_signal():
from scipy import signal
import numpy as np
b, a = signal.butter(4, 0.1)
filtered = signal.filtfilt(b, a, np.random.randn(1000))
return f"フィルタ出力 len={len(filtered)}"
def _test_interpolate():
from scipy.interpolate import CubicSpline
import numpy as np
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 0, 1, 0])
cs = CubicSpline(x, y)
return f"補間値 f(2.5)={cs(2.5):.3f}"
def _test_linalg():
from scipy import linalg
import numpy as np
A = np.random.rand(100, 100)
vals = linalg.eigvals(A)
return f"固有値 {len(vals)}個"
def _test_integrate():
from scipy import integrate
import numpy as np
result, _ = integrate.quad(lambda x: np.exp(-x**2), -np.inf, np.inf)
return f"∫e^(-x²) = {result:.4f}"
def _test_sparse():
from scipy import sparse
m = sparse.random(1000, 1000, density=0.01, format="csr")
return f"疎行列 nnz={m.nnz}"
if __name__ == "__main__":
diagnose_scipy()
11. ユースケース別ガイド
ユースケース1: 実験データの曲線フィッティングとパラメータ推定
from scipy.optimize import curve_fit
from scipy import stats
import numpy as np
# 酵素反応速度のミカエリス・メンテン式: V = Vmax * [S] / (Km + [S])
def michaelis_menten(S, Vmax, Km):
return Vmax * S / (Km + S)
# 実験データ（基質濃度と反応速度）
S_data = np.array([0.5, 1, 2, 5, 10, 20, 50, 100])
V_data = np.array([2.1, 3.8, 5.9, 8.5, 9.8, 10.8, 11.3, 11.6])
V_err = np.array([0.3, 0.3, 0.4, 0.3, 0.4, 0.5, 0.3, 0.4]) # 測定誤差
# 重み付きカーブフィッティング（誤差が小さいデータ点ほど重視）
popt, pcov = curve_fit(michaelis_menten, S_data, V_data, sigma=V_err, p0=[12, 3])
Vmax_fit, Km_fit = popt
perr = np.sqrt(np.diag(pcov))
print(f"Vmax = {Vmax_fit:.2f} ± {perr[0]:.2f}")
print(f"Km = {Km_fit:.2f} ± {perr[1]:.2f}")
# フィッティングの良さを評価（残差の正規性検定）
residuals = V_data - michaelis_menten(S_data, *popt)
_, p_normal = stats.shapiro(residuals)
print(f"残差の正規性検定 p値: {p_normal:.4f}")
ユースケース2: 時系列データのノイズ除去とピーク検出
from scipy import signal
import numpy as np
# 心拍数データのシミュレーション（60BPM = 1Hz）
fs = 250 # サンプリング周波数
t = np.arange(0, 10, 1/fs)
ecg_clean = np.sin(2 * np.pi * 1.0 * t) # 1Hzの心拍
noise = 0.3 * np.random.randn(len(t)) + 0.1 * np.sin(2 * np.pi * 50 * t) # ノイズ+電源ハム
ecg_noisy = ecg_clean + noise
# バンドパスフィルタ（0.5〜30Hzを通過）
sos = signal.butter(4, [0.5, 30], btype="bandpass", fs=fs, output="sos")
ecg_filtered = signal.sosfiltfilt(sos, ecg_noisy)
# ピーク検出（心拍のR波）
peaks, properties = signal.find_peaks(ecg_filtered, distance=fs*0.6, height=0.5)
heart_rate = 60 * fs / np.diff(peaks).mean() if len(peaks) > 1 else 0
print(f"検出されたピーク数: {len(peaks)}")
print(f"推定心拍数: {heart_rate:.0f} BPM")
ユースケース3: 大規模疎行列でのページランク計算
from scipy import sparse
from scipy.sparse.linalg import eigs
import numpy as np
# Webページのリンク構造を疎行列で表現
np.random.seed(42)
n_pages = 10000
n_links = 50000 # 各ページから平均5リンク
rows = np.random.randint(0, n_pages, n_links)
cols = np.random.randint(0, n_pages, n_links)
data = np.ones(n_links)
# 隣接行列（疎行列）
A = sparse.csr_matrix((data, (rows, cols)), shape=(n_pages, n_pages))
# 列の正規化（確率遷移行列に変換）
col_sums = np.array(A.sum(axis=0)).flatten()
col_sums[col_sums == 0] = 1 # ゼロ除算防止
D_inv = sparse.diags(1.0 / col_sums)
M = A @ D_inv
# べき乗法でPageRankベクトルを計算
rank = np.ones(n_pages) / n_pages
damping = 0.85
for _ in range(100):
rank = damping * M @ rank + (1 - damping) / n_pages
# 上位ページ
top_pages = np.argsort(rank)[::-1][:5]
print("PageRank 上位5ページ:")
for i, page in enumerate(top_pages):
print(f" {i+1}位: ページ{page} (スコア: {rank[page]:.6f})")
print(f"\n疎行列メモリ: {(A.data.nbytes + A.indices.nbytes + A.indptr.nbytes) / 1024**2:.1f} MB")
print(f"密行列なら: {n_pages**2 * 8 / 1024**2:.0f} MB")
12. 学習ロードマップ
| レベル | 目安期間 | 到達状態 |
|---|---|---|
| Lv.1 基礎 | 1〜2週間 | 最適化と統計検定の基本が使える |
| Lv.2 実践 | 2〜4週間 | 信号処理、補間、線形代数を実務データに適用できる |
| Lv.3 応用 | 1〜2ヶ月 | 微分方程式の数値解法、疎行列を使った大規模計算ができる |
| Lv.4 発展 | 2ヶ月〜 | 記号計算やML/DLとの連携で研究・開発レベルの計算ができる |
まとめ
この記事では、SciPyの全体像を「大学の理工学部棟」の比喩で解き明かしてきた。
SciPyはNumPyという共通インフラの上に建つ専門研究室の集合体だ。最適化研究室はコスト関数を最小化し、統計学研究室はデータの差が偶然か否かを判定し、信号処理研究室はノイズまみれのデータから本来の信号を取り出す。補間研究室は離散データの隙間を埋め、線形代数研究室は大規模な連立方程式を解き、積分研究室は解析解が存在しない積分を数値的に計算する。そして疎行列研究室は、ほとんどゼロの巨大行列を99%のメモリ削減で扱う魔法を提供する。
個人的な体験として、筆者がローカルLLMの応答品質を統計的に評価しようとした際、「モデルAとモデルBの回答品質スコアに有意差があるか」を scipy.stats.mannwhitneyu で検定した。正規分布を仮定できないスコアデータだったので、ノンパラメトリック検定が必要だった。SciPyがなければ自前で実装するか、R言語に逃げるかの二択だったと思う...orz
SciPyは「NumPyの上位互換」ではなく、NumPyでは手が届かない専門領域への扉だ。最適化、統計、信号処理、補間、線形代数、積分、疎行列。これらの研究室がすべて同じ建物内にあり、NumPy配列という共通言語でデータをやり取りできる。この統一性こそが、SciPyが科学技術計算のデファクトスタンダードであり続ける理由だ。
参考文献
Python・NumPy・pandasの基礎から学びたい方は、こちらのシリーズもどうぞ:
筆者のXアカウントはこちら（@geneLab_999）