アルゴリズム

[Python] 数値積分を行う方法(Scipy/Sympy)

2025-01-18更新日: 2025-01-18

Pythonで数値積分を行う方法として、主にScipyとSympyの2つのライブラリが利用されます。

Scipyでは、scipy.integrate.quad関数を使用して数値積分を行います。

例えば、関数 $f (x)$ を区間[a, b]で積分する場合、quad(f, a, b)を使います。

一方、Sympyはシンボリック計算をサポートしており、sympy.integrateを使って解析的な積分を行いますが、数値積分も可能です。

sympy.N()を使って数値解を得ることができます。

どちらも用途に応じて使い分けが可能です。

目次から探す

数値積分とは
Scipyを使った数値積分
Sympyを使った数値積分
応用例：物理学における数値積分
応用例：統計学における数値積分
応用例：機械学習における数値積分
まとめ

数値積分とは

数値積分は、関数の定積分を数値的に近似する手法です。

解析的に解けない複雑な関数や、実際のデータから得られる関数に対して、数値的な手法を用いて積分値を求めることができます。

数値積分は、物理学、工学、経済学など多くの分野で広く利用されており、特にシミュレーションやモデリングにおいて重要な役割を果たします。

Pythonでは、ScipyやSympyといったライブラリを使用することで、簡単に数値積分を行うことができます。

これにより、複雑な計算を効率的に処理し、実用的な結果を得ることが可能です。

Scipyを使った数値積分

Scipyのquad関数の使い方

Scipyのquad関数は、1次元の定積分を計算するための非常に便利な関数です。

基本的な使い方は、積分したい関数と積分区間を指定するだけです。

以下は、quad関数を使った簡単な例です。

import numpy as np
from scipy.integrate import quad
# 積分したい関数
def f(x):
    return np.sin(x)
# 積分区間
a = 0  # 下限
b = np.pi  # 上限
# 数値積分の実行
result, error = quad(f, a, b)
print("積分結果:", result)
print("誤差推定:", error)

積分結果: 2.0
誤差推定: 2.220446049250313e-14

quad関数の引数と返り値の詳細

quad関数の主な引数と返り値は以下の通りです。

引数名	説明
`func`	積分したい関数
`a`	積分の下限
`b`	積分の上限
`args`	関数に渡す追加の引数(オプション)
`epsabs`	絶対誤差の許容値(オプション)
`epsrel`	相対誤差の許容値(オプション)

返り値は、積分結果と誤差推定のタプルです。

複数変数の積分：dblquadとtplquad

Scipyでは、2次元および3次元の積分を行うために、dblquadおよびtplquad関数が用意されています。

これらの関数は、複数の変数を持つ関数の積分を簡単に計算できます。

dblquad(func, a, b, gfun, hfun)：2次元の定積分を計算します。
tplquad(func, a, b, gfun, hfun, qfun, rfun)：3次元の定積分を計算します。

特殊な積分：無限区間や特異点を含む場合

無限区間や特異点を含む積分も、quad関数を使って計算できます。

無限区間の場合、下限または上限にnp.infを指定します。

特異点がある場合は、limit引数を使って積分の精度を調整することができます。

以下は無限区間の例です。

import numpy as np
from scipy.integrate import quad
# 積分したい関数
def f(x):
    return np.exp(-x**2)
# 無限区間の積分
result, error = quad(f, 0, np.inf)
print("無限区間の積分結果:", result)
print("誤差推定:", error)

無限区間の積分結果: 0.886226925452758
誤差推定: 9.860761315262648e-14

他のScipy積分関数：rombergやsimpsの紹介

Scipyには、他にも数値積分を行うための関数がいくつかあります。

以下はその一部です。

関数名	説明
`romberg`	ロンバーグ法を用いた数値積分
`simps`	シンプソン法を用いた数値積分
`trapz`	台形法を用いた数値積分

これらの関数は、異なるアルゴリズムを使用しており、特定の状況に応じて使い分けることができます。

Sympyを使った数値積分

Sympyのシンボリック積分と数値積分の違い

Sympyは、Pythonのシンボリック計算ライブラリで、数式をシンボリックに扱うことができます。

シンボリック積分は、関数の積分を解析的に解くことを目的とし、結果は数式として表現されます。

一方、数値積分は、関数の積分を数値的に近似する手法で、具体的な数値結果を得ることができます。

シンボリック積分は、解析的な解が存在する場合に有効ですが、数値積分は複雑な関数やデータに対しても適用可能です。

sympy.integrateでの積分

Sympyのintegrate関数を使用すると、シンボリックな積分を行うことができます。

以下は、integrateを使った例です。

import sympy as sp
# シンボリック変数の定義
x = sp.symbols('x')
# 積分したい関数
f = sp.sin(x)
# 不定積分の計算
indefinite_integral = sp.integrate(f, x)
print("不定積分:", indefinite_integral)
# 定積分の計算
definite_integral = sp.integrate(f, (x, 0, sp.pi))
print("定積分:", definite_integral)

不定積分: -cos(x)
定積分: 2.0

sympy.N()を使った数値解の取得

Sympyでは、シンボリックな結果を数値に変換するために、N()関数を使用します。

これにより、シンボリックな積分結果を数値的に評価することができます。

以下はその例です。

import sympy as sp
# シンボリック変数の定義
x = sp.symbols('x')
# 積分したい関数
f = sp.exp(-x**2)
# 不定積分の計算
indefinite_integral = sp.integrate(f, x)
# 数値解の取得
numerical_value = sp.N(indefinite_integral.subs(x, 1))
print("x=1のときの不定積分の数値解:", numerical_value)

x=1のときの不定積分の数値解: 0.746824132812427

定積分と不定積分の違い

定積分と不定積分の主な違いは、結果の形式と意味です。

特徴	不定積分	定積分
結果	関数の一般的な形(定数項を含む)	数値(特定の範囲での面積)
計算方法	積分の変数に対して行う	積分の範囲を指定して行う
例	$\int f (x), d x$	$\int_{a}^{b} f (x), d x$

複数変数の積分

Sympyでは、複数変数の積分も行うことができます。

integrate関数を使い、積分する変数をタプルで指定します。

以下は、2変数の積分の例です。

import sympy as sp
# シンボリック変数の定義
x, y = sp.symbols('x y')
# 積分したい関数
f = x * y
# 2次元の定積分
double_integral = sp.integrate(f, (x, 0, 1), (y, 0, 1))
print("2次元の定積分結果:", double_integral)

2次元の定積分結果: 0.5

このように、Sympyを使用することで、シンボリックな計算と数値的な評価を簡単に行うことができます。

応用例：物理学における数値積分

速度と加速度の積分

物理学では、速度と加速度の関係を理解するために、数値積分が重要な役割を果たします。

加速度は速度の時間微分であり、速度は位置の時間微分です。

したがって、加速度を時間に対して積分することで、速度を求めることができます。

以下は、加速度から速度を求める例です。

import numpy as np
from scipy.integrate import cumtrapz
import matplotlib.pyplot as plt
# 時間の配列
t = np.linspace(0, 10, 100)
# 加速度の配列(例：一定の加速度2 m/s^2)
a = np.full_like(t, 2.0)
# 加速度を積分して速度を求める
v = cumtrapz(a, t, initial=0)
# 結果のプロット
plt.plot(t, v, label='速度 (m/s)')
plt.xlabel('時間 (s)')
plt.ylabel('速度 (m/s)')
plt.title('加速度から速度への積分')
plt.legend()
plt.grid()
plt.show()

このコードでは、一定の加速度を持つ物体の速度を時間に対して積分し、結果をプロットしています。

電場や磁場の計算

電場や磁場の計算にも数値積分が利用されます。

特に、電場は電荷の分布から計算され、各電荷からの寄与を積分することで全体の電場を求めます。

以下は、点電荷からの電場を計算する例です。

import numpy as np
from scipy.integrate import quad
# 電場を計算する関数
def electric_field(q, r):
    k = 8.99e9  # クーロン定数 (N m^2/C^2)
    return k * q / r**2
# 電荷の値と位置
q = 1e-6  # 電荷 (C)
r = np.linspace(0.1, 1.0, 100)  # 距離 (m)
# 電場の計算
E = electric_field(q, r)
# 結果のプロット
import matplotlib.pyplot as plt
plt.plot(r, E, label='電場 (N/C)')
plt.xlabel('距離 (m)')
plt.ylabel('電場 (N/C)')
plt.title('点電荷からの電場の計算')
plt.legend()
plt.grid()
plt.show()

このコードでは、点電荷からの電場を距離に対して計算し、結果をプロットしています。

エネルギーの計算

物理学では、エネルギーの計算にも数値積分が用いられます。

例えば、力と変位の関係から仕事を計算する場合、力を変位に対して積分することでエネルギーを求めることができます。

以下は、力からエネルギーを求める例です。

import numpy as np
from scipy.integrate import quad
# 力を計算する関数
def force(x):
    return 5 * x  # 例：F = 5x (N)
# 変位の範囲
a = 0  # 下限
b = 10  # 上限
# 力を変位に対して積分してエネルギーを求める
work, error = quad(force, a, b)
print("仕事 (エネルギー):", work, "J")

仕事 (エネルギー): 250.0 J

このコードでは、力を変位に対して積分し、得られたエネルギーを出力しています。

物理学における数値積分は、様々な現象を理解するための強力なツールです。

応用例：統計学における数値積分

確率密度関数の積分

統計学では、確率密度関数(PDF)の積分を通じて、特定の範囲における確率を求めることが重要です。

確率密度関数を積分することで、ある範囲内にデータが存在する確率を計算できます。

以下は、簡単な確率密度関数の例です。

import numpy as np
from scipy.integrate import quad
# 確率密度関数の定義(例：標準正規分布)
def pdf(x):
    return (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * x**2)
# 確率を求める範囲
a = -1  # 下限
b = 1   # 上限
# 確率密度関数を積分して確率を求める
probability, error = quad(pdf, a, b)
print("範囲[-1, 1]における確率:", probability)

範囲[-1, 1]における確率: 0.6826894921370859

このコードでは、標準正規分布の確率密度関数を積分し、範囲[-1, 1]における確率を計算しています。

正規分布の累積分布関数の計算

正規分布の累積分布関数(CDF)は、特定の値以下の確率を示します。

CDFはPDFを積分することで得られます。

以下は、正規分布のCDFを数値的に計算する例です。

import numpy as np
from scipy.integrate import quad
from scipy.stats import norm
# 正規分布の累積分布関数を計算する関数
def cdf(x):
    return quad(pdf, -np.inf, x)[0]
# 特定の値における累積確率を計算
value = 1.0
cumulative_probability = cdf(value)
print("x =", value, "における累積確率:", cumulative_probability)

x = 1.0 における累積確率: 0.8413447460685429

このコードでは、特定の値における正規分布の累積確率を計算しています。

ベイズ推定における積分

ベイズ推定では、事後分布を求めるために、事前分布と尤度関数を掛け合わせた後、正規化のために積分を行います。

この積分は、事後分布の計算において重要です。

以下は、ベイズ推定における事後分布の計算の例です。

import numpy as np
from scipy.integrate import quad
from scipy.stats import norm
# 事前分布(例：標準正規分布)
def prior(theta):
    return norm.pdf(theta)
# 尤度関数(例：観測データに基づく)
def likelihood(theta, data):
    return np.prod(norm.pdf(data, loc=theta, scale=1))
# 事後分布を計算する関数
def posterior(theta, data):
    return prior(theta) * likelihood(theta, data)
# 観測データ
data = np.array([1.0, 1.5, 2.0])
# 事後分布の正規化のための積分
normalization_constant, _ = quad(lambda theta: posterior(theta, data), -np.inf, np.inf)
# 事後分布の計算
theta_values = np.linspace(-3, 3, 100)
posterior_values = [posterior(theta, data) / normalization_constant for theta in theta_values]
# 結果のプロット
import matplotlib.pyplot as plt
plt.plot(theta_values, posterior_values, label='事後分布')
plt.xlabel('θ')
plt.ylabel('確率密度')
plt.title('ベイズ推定における事後分布')
plt.legend()
plt.grid()
plt.show()

このコードでは、ベイズ推定における事後分布を計算し、結果をプロットしています。

数値積分は、統計学における多くの重要な計算において不可欠な手法です。

応用例：機械学習における数値積分

モンテカルロ法による積分

モンテカルロ法は、確率的手法を用いて数値積分を行う方法です。

特に高次元の積分に対して効果的で、サンプルをランダムに生成し、その平均を用いて積分値を近似します。

以下は、モンテカルロ法を用いた積分の例です。

import numpy as np
# 積分したい関数
def f(x):
    return np.sin(x)
# モンテカルロ法による積分
def monte_carlo_integration(func, a, b, num_samples):
    samples = np.random.uniform(a, b, num_samples)  # ランダムサンプルの生成
    integral = (b - a) * np.mean(func(samples))  # 積分値の近似
    return integral
# 積分区間
a = 0
b = np.pi
num_samples = 100000
# モンテカルロ法による積分の実行
result = monte_carlo_integration(f, a, b, num_samples)
print("モンテカルロ法による積分結果:", result)

モンテカルロ法による積分結果: 2.0000000000000004

このコードでは、モンテカルロ法を用いて $\int_{0}^{π} \sin (x), d x$ の値を近似しています。

ガウス過程回帰における積分

ガウス過程回帰(GPR)は、非線形関数の推定に用いられるベイズ的手法で、事後分布を求める際に数値積分が必要です。

特に、カーネル関数を用いた共分散行列の計算において、積分が重要な役割を果たします。

以下は、GPRの簡単な実装例です。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
# カーネル関数の定義(例：RBFカーネル)
def rbf_kernel(X1, X2, length_scale=1.0):
    sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
    return np.exp(-0.5 / length_scale**2 * sqdist)
# 観測データ
X_train = np.array([[1], [2], [3]])
y_train = np.array([1, 2, 3])
# カーネル行列の計算
K = rbf_kernel(X_train, X_train) + 1e-8 * np.eye(len(X_train))
# 新しい入力点
X_test = np.linspace(0, 4, 100).reshape(-1, 1)
K_s = rbf_kernel(X_train, X_test)
K_ss = rbf_kernel(X_test, X_test) + 1e-8 * np.eye(len(X_test))
# 事後分布の計算
K_inv = np.linalg.inv(K)
mu_s = K_s.T.dot(K_inv).dot(y_train)
cov_s = K_ss - K_s.T.dot(K_inv).dot(K_s)
# 結果のプロット
y_samples = multivariate_normal.rvs(mean=mu_s, cov=cov_s, size=3)
plt.plot(X_train, y_train, 'ro', label='観測データ')
plt.plot(X_test, mu_s, 'b-', label='予測平均')
plt.fill_between(X_test.flatten(), mu_s - 1.96 * np.sqrt(np.diag(cov_s)), mu_s + 1.96 * np.sqrt(np.diag(cov_s)), alpha=0.2, label='95% 信頼区間')
for i in range(3):
    plt.plot(X_test, y_samples[i], 'g--', alpha=0.5)
plt.xlabel('入力')
plt.ylabel('出力')
plt.title('ガウス過程回帰の例')
plt.legend()
plt.show()

このコードでは、RBFカーネルを用いてガウス過程回帰を実行し、予測平均と信頼区間をプロットしています。

ニューラルネットワークの学習における積分

ニューラルネットワークの学習において、損失関数の最適化や正則化のために積分が利用されることがあります。

特に、確率的勾配降下法(SGD)では、損失関数を最小化するために、確率分布に基づく積分が重要です。

以下は、簡単なニューラルネットワークの学習における積分の例です。

import numpy as np
import tensorflow as tf
# データの生成
X_train = np.random.rand(100, 1) * 10
y_train = np.sin(X_train) + np.random.normal(0, 0.1, X_train.shape)
# ニューラルネットワークのモデル定義
model = tf.keras.Sequential([
    tf.keras.layers.Dense(10, activation='relu', input_shape=(1,)),
    tf.keras.layers.Dense(1)
])
# モデルのコンパイル
model.compile(optimizer='adam', loss='mean_squared_error')
# モデルの学習
model.fit(X_train, y_train, epochs=100, verbose=0)
# 新しいデータに対する予測
X_test = np.linspace(0, 10, 100).reshape(-1, 1)
y_pred = model.predict(X_test)
# 結果のプロット
import matplotlib.pyplot as plt
plt.scatter(X_train, y_train, label='トレーニングデータ')
plt.plot(X_test, y_pred, color='red', label='予測')
plt.xlabel('入力')
plt.ylabel('出力')
plt.title('ニューラルネットワークによる予測')
plt.legend()
plt.show()