SciPyを使ったPythonでの科学計算

SciPyとは

SciPyは、科学技術計算を効率的に行うためのPythonライブラリで、NumPyを基盤として構築されています。
数値計算や統計、信号処理、線形代数、最適化など、様々なタスクを効率的に処理することができます。

インストール

SciPyはPythonの外部ライブラリなので、まずはターミナルでライブラリをインストールする必要があります。

pip install scipy

SciPyのほぼすべての機能がNumPyのndarray（多次元配列）に依存しているので、SciPyをインストールすると、同時にNumPyも自動的にインストールされます。

基本的な使い方

SciPyは大きく以下のようなサブパッケージに分かれています。それぞれのサブパッケージには、特定の計算分野に対応した関数やクラスが集められています。

サブパッケージ名	説明
scipy.linalg	線形代数の関数（行列計算、行列の逆行列、固有値など）
scipy.optimize	最適化と方程式の解法
scipy.stats	統計的な分布、仮説検定、確率計算
scipy.integrate	数値積分と常微分方程式の解法
scipy.signal	信号処理（フィルタリング、変換、スペクトル解析など）
scipy.fftpack	高速フーリエ変換と逆変換
scipy.spatial	空間データの計算（距離計算、ボロノイ図など）
scipy.interpolate	データの補間（スプライン補間、多次元補間）
scipy.ndimage	TN次元の画像処理his

これらのサブパッケージを適切に使い分けることで、さまざまな数値計算や科学技術計算を効率的に行うことができます。

1.線形代数（scipy.linalg）

SciPyの線形代数サブパッケージは、行列計算、逆行列の計算、固有値問題の解法など、さまざまな行列演算を効率的に行うことができます。

・逆行列の計算

行列の逆行列を計算するには、linalg.inv()関数を使用します。

import numpy as np
from scipy import linalg

# 2x2の行列を定義
A = np.array([[1, 2],
              [3, 4]])

# 逆行列を計算
inv_A = linalg.inv(A)

print("逆行列：\n", inv_A)

# 出力
逆行列:
 [[-2.   1. ]
 [ 1.5 -0.5]]

このコードでは、linalg.inv()を使用して行列Aの逆行列を計算しています。行列が正則でない場合（つまり、逆行列が存在しない場合）、関数はエラーを発生させます。

・行列の固有値と固有ベクトル

SciPyでは、行列の固有値と固有ベクトルも簡単に計算できます。これにはlinalg.eig()を使用します。

# 固有値と固有ベクトルを計算
eigenvalues, eigenvectors = linalg.eig(A)

print("固有値：\n", eigenvalues)
print("固有ベクトル：\n", eigenvectors)

# 出力
固有値:
 [-0.37228132+0.j  5.37228132+0.j]
固有ベクトル:
 [[-0.82456484 -0.41597356]
  [ 0.56576746 -0.90937671]]

2.最適化（scipy.optimize）

最適化は、ある関数の最小値や最大値を見つけることを目的とした計算です。scipy.optimizeには、関数の最小化や制約付きの最適化、方程式の解を見つけるためのアルゴリズムが含まれています。

・目的関数の最小化

optimize.minimize()を使って、関数の最小値を求めることができます。

from scipy import optimize
import numpy as np

# 目的関数を定義
def f(x):
    return x**2 + 5 * np.sin(x)

# 初期値x0=0で最小化を開始
result = optimize.minimize(f, x0=0)

print("最小値を取るxの値：", result.x)
print("目的関数の最小値：", result.fun)

# 出力
最小値を取るxの値: [-1.11051058]
目的関数の最小値: -3.2463942726911937

ここでは、f(x)という関数（x**2 + 5 * np.sin(x)）を定義し、その最小値を見つけるためにminimize()を使っています。初期値としてx=0から探索を開始し、最小値が見つかる場所とその値を出力しています。

方程式の解を見つける

非線形方程式の解を見つけるには、optimize.root()を使用します。

# 方程式 cos(x) = x の解を探す
def equation(x):
    return np.cos(x) - x

solution = optimize.root(equation, x0=0)

print("方程式の解：", solution.x)

# 出力
方程式の解: [0.73908513]

ここでは、np.cos(x) - x = 0という方程式を解くためにoptimize.root()を使っています。x0=0は初期推定値で、cos(x) = xの解としてx ≈ 0.739が見つかります。

3.統計解析（scipy.stats）

統計解析を行うために、SciPyのstatsサブパッケージを使用します。これにより、統計的な検定、確率分布、サンプルデータの生成などが簡単に行えます。

・正規分布のサンプル生成

SciPyを使って、正規分布に従う乱数データを生成し、その基本統計量を計算する例です。

from scipy import stats

# 平均0、標準偏差1の正規分布に従うサンプルを1000個生成
data = stats.norm.rvs(size=1000, loc=0, scale=1)

# 平均と標準偏差を計算
mean = np.mean(data)
std_dev = np.std(data)

print("データの平均：", mean)
print("データの標準偏差：", std_dev)

# 出力
データの平均: -0.03278572502910955
データの標準偏差: 0.9821633212835997

・t検定

SciPyを使って、1つのサンプルデータに対するt検定を行うことができます。これにより、サンプルが特定の平均値に従うかどうかを検定します。

# t検定を実施（平均0かどうか）
t_stat, p_value = stats.ttest_1samp(data, 0)

print("t検定のt値：", t_stat)
print("p値：", p_value)

# 出力
t検定のt値: -1.0527051531751178
p値: 0.29277281056901965

このt検定では、データが平均0の正規分布に従うかどうかを検定しています。p値が0.05以上であるため、このサンプルは平均0に従うという仮説を棄却できません。

4.信号処理（scipy.signal）

SciPyのsignalサブパッケージは、信号処理に特化した機能を提供しています。フィルタリングやフーリエ変換、信号スペクトルの解析に使われます。

・信号のフィルタリング

SciPyを使って信号にローパスフィルタを適用し、高周波成分を除去します。

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

# サンプル信号を作成
t = np.linspace(0, 1.0, 200)
x = np.sin(2 * np.pi * 5 * t) + np.sin(2 * np.pi * 50 * t)

# ローパスフィルタを設計
b, a = signal.butter(4, 0.1)
y = signal.filtfilt(b, a, x)

# 信号のプロット
plt.plot(t, x, label='Noisy signal')
plt.plot(t, y, label='Filtered signal')
plt.legend()
plt.show()

このコードでは、50Hzのノイズが含まれた信号を生成し、ローパスフィルタを適用してノイズを除去しています。フィルタリング前後の信号がプロットされています。

応用例

微分方程式の解法（scipy.integrate）

SciPyのintegrateサブパッケージでは、数値積分や微分方程式の解法が簡単に行えます。

常微分方程式の解法

シンプルな常微分方程式（dy/dt = -y）をSciPyで数値的に解く例を示します。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# 微分方程式 dy/dt = -y を定義
def model(t, y):
    return -y

# t=0からt=5までの間を解く
solution = solve_ivp(model, [0, 5], [1], t_eval=np.linspace(0, 5, 100))

# 解のプロット
plt.plot(solution.t, solution.y[0])
plt.xlabel('Time')
plt.ylabel('y')
plt.show()

初期値y(0) = 1のシンプルな指数関数的減衰を表す微分方程式を解いています。solve_ivpは、数値的に微分方程式を解くための強力な関数で、さまざまな方程式に対応しています。