Pythonは「初心者にやさしい」だけでなく、「数学処理にも強い」言語です。
簡単な計算から、三角関数、統計、行列計算、微分・積分まで幅広く対応でき、教育・研究・実務の現場でも大活躍しています。
Pythonが数学処理に強い理由
- 豊富な標準ライブラリ
- NumPy、SciPyなどの高性能な外部ライブラリ
- 直感的でわかりやすい文法
- 科学計算コミュニティが活発
この記事では、Pythonでよく使われる数学処理の方法を体系的に整理し、標準ライブラリや代表的な外部ライブラリの使い方も合わせて紹介します。
基本の数値計算（演算子と標準関数）
四則演算と基本演算子
Pythonでは、様々な数値演算を直感的に行えます。
# 基本的な演算
print(3 + 2) # 足し算: 5
print(5 - 1) # 引き算: 4
print(4 * 2) # 掛け算: 8
print(6 / 2) # 割り算: 3.0（浮動小数点）
print(7 // 2) # 切り捨て割り算: 3
print(7 % 2) # 余り: 1
print(2 ** 3) # 累乗: 8
比較演算子
# 比較演算
print(5 > 3) # より大きい: True
print(5 < 3) # より小さい: False
print(5 >= 5) # 以上: True
print(5 <= 4) # 以下: False
print(5 == 5) # 等しい: True
print(5 != 3) # 等しくない: True
数値型とその特徴
# 整数型（int）
a = 10
print(type(a)) # <class 'int'>
# 浮動小数点型（float）
b = 3.14
print(type(b)) # <class 'float'>
# 複素数型（complex）
c = 3 + 4j
print(type(c)) # <class 'complex'>
print(c.real) # 実部: 3.0
print(c.imag) # 虚部: 4.0
組み込み数学関数
# 絶対値
print(abs(-5)) # 5
print(abs(3 + 4j)) # 複素数の絶対値: 5.0
# 最大値・最小値
print(max(1, 5, 3)) # 5
print(min([2, 8, 1])) # 1
# 合計
print(sum([1, 2, 3, 4])) # 10
# 累乗
print(pow(2, 3)) # 8
print(pow(2, 3, 5)) # (2^3) % 5 = 3
# 四捨五入
print(round(3.14159, 2)) # 3.14
# 商と余り
print(divmod(17, 5)) # (3, 2)
標準ライブラリ math の活用
mathモジュールは、数学的な関数と定数を提供する標準ライブラリです。
import math
基本的な数学関数
# 平方根
print(math.sqrt(16)) # 4.0
print(math.sqrt(2)) # 1.4142135623730951
# 累乗と指数
print(math.pow(2, 3)) # 8.0
print(math.exp(1)) # e^1 = 2.718281828459045
# 対数
print(math.log(10)) # 自然対数: 2.302585092994046
print(math.log10(100)) # 常用対数: 2.0
print(math.log(8, 2)) # 底2の対数: 3.0
print(math.log2(8)) # 底2の対数: 3.0
# 階乗
print(math.factorial(5)) # 5! = 120
三角関数
# 基本的な三角関数（ラジアン）
print(math.sin(math.pi / 2)) # sin(π/2) = 1.0
print(math.cos(0)) # cos(0) = 1.0
print(math.tan(math.pi / 4)) # tan(π/4) = 1.0
# 逆三角関数
print(math.asin(1)) # arcsin(1) = π/2
print(math.acos(0)) # arccos(0) = π/2
print(math.atan(1)) # arctan(1) = π/4
print(math.atan2(1, 1)) # atan2(y, x) = π/4
# 双曲線関数
print(math.sinh(0)) # sinh(0) = 0.0
print(math.cosh(0)) # cosh(0) = 1.0
print(math.tanh(0)) # tanh(0) = 0.0
角度変換
# 度とラジアンの変換
print(math.radians(180)) # 180度 → π ラジアン
print(math.degrees(math.pi)) # π ラジアン → 180度
# 実用例：度で計算
angle_deg = 30
angle_rad = math.radians(angle_deg)
print(f"sin({angle_deg}°) = {math.sin(angle_rad)}") # sin(30°) = 0.5
数学定数
# 重要な数学定数
print(f"円周率 π = {math.pi}") # 3.141592653589793
print(f"ネイピア数 e = {math.e}") # 2.718281828459045
print(f"黄金比 φ = {(1 + math.sqrt(5)) / 2}") # 1.618033988749895
# 無限大とNaN
print(math.inf) # inf（無限大）
print(math.nan) # nan（非数）
print(math.isnan(math.nan)) # True
print(math.isinf(math.inf)) # True
床関数・天井関数
# 床関数（切り下げ）
print(math.floor(3.7)) # 3
print(math.floor(-3.7)) # -4
# 天井関数（切り上げ）
print(math.ceil(3.2)) # 4
print(math.ceil(-3.2)) # -3
# 整数部分の切り出し
print(math.trunc(3.7)) # 3
print(math.trunc(-3.7)) # -3
最大公約数・最小公倍数
# 最大公約数
print(math.gcd(48, 18)) # 6
print(math.gcd(24, 36, 48)) # 12（Python 3.9以降）
# 最小公倍数（Python 3.9以降）
print(math.lcm(12, 18)) # 36
print(math.lcm(4, 6, 8)) # 24
統計処理（statisticsモジュール）
statisticsモジュールは、基本的な統計計算を提供します。
import statistics as stat
中心傾向の測度
data = [10, 20, 30, 40, 50]
# 平均値
print(stat.mean(data)) # 30.0
# 調和平均
print(stat.harmonic_mean(data)) # 22.136515775282186
# 幾何平均
print(stat.geometric_mean(data)) # 25.992104989487313
# 中央値
print(stat.median(data)) # 30
# 最頻値
data_with_mode = [1, 2, 2, 3, 4]
print(stat.mode(data_with_mode)) # 2
# 複数の最頻値
data_multimode = [1, 1, 2, 2, 3]
print(stat.multimode(data_multimode)) # [1, 2]
散布度の測度
data = [2, 4, 6, 8, 10]
# 分散
print(stat.variance(data)) # 10.0（標本分散）
print(stat.pvariance(data)) # 8.0（母分散）
# 標準偏差
print(stat.stdev(data)) # 3.1622776601683795（標本標準偏差）
print(stat.pstdev(data)) # 2.8284271247461903（母標準偏差）
# 分位数
data_large = list(range(1, 101)) # 1から100まで
print(stat.quantiles(data_large, n=4)) # 四分位数
実用的な統計例
# 成績データの分析
scores = [85, 90, 78, 92, 88, 76, 95, 82, 89, 91]
print(f"平均点: {stat.mean(scores):.1f}")
print(f"中央値: {stat.median(scores)}")
print(f"標準偏差: {stat.stdev(scores):.2f}")
print(f"最高点: {max(scores)}")
print(f"最低点: {min(scores)}")
# 偏差値の計算
def calculate_deviation_score(score, scores):
mean_val = stat.mean(scores)
std_val = stat.stdev(scores)
return 50 + (score - mean_val) / std_val * 10
for score in scores:
deviation = calculate_deviation_score(score, scores)
print(f"点数: {score}, 偏差値: {deviation:.1f}")
NumPyによる高速な数学処理
NumPyは科学計算用ライブラリで、ベクトル・行列・配列計算が超高速に行えます。
import numpy as np
配列の作成と基本操作
# 配列の作成
a = np.array([1, 2, 3, 4, 5])
b = np.array([[1, 2], [3, 4]])
print(a) # [1 2 3 4 5]
print(b) # [[1 2]
# [3 4]]
# 特殊な配列
zeros = np.zeros(5) # [0. 0. 0. 0. 0.]
ones = np.ones((2, 3)) # 2×3の1で埋められた配列
eye = np.eye(3) # 3×3の単位行列
arange_arr = np.arange(0, 10, 2) # [0 2 4 6 8]
linspace_arr = np.linspace(0, 1, 5) # [0. 0.25 0.5 0.75 1. ]
ベクトル演算
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
# 要素ごとの演算
print(a + b) # [5 7 9]
print(a - b) # [-3 -3 -3]
print(a * b) # [4 10 18]
print(a / b) # [0.25 0.4 0.5 ]
# ベクトルの内積
print(np.dot(a, b)) # 32
print(a @ b) # 32（Python 3.5以降）
# ベクトルの外積（3次元）
print(np.cross(a, b)) # [-3 6 -3]
# ベクトルのノルム
print(np.linalg.norm(a)) # 3.7416573867739413
行列演算
# 行列の作成
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
# 行列の加減算
print(A + B) # [[ 6 8]
# [10 12]]
# 行列の積
print(A @ B) # [[19 22]
# [43 50]]
# 転置行列
print(A.T) # [[1 3]
# [2 4]]
# 逆行列
print(np.linalg.inv(A)) # [[-2. 1. ]
# [ 1.5 -0.5]]
# 行列式
print(np.linalg.det(A)) # -2.0000000000000004
# 固有値と固有ベクトル
eigenvalues, eigenvectors = np.linalg.eig(A)
print("固有値:", eigenvalues)
print("固有ベクトル:", eigenvectors)
数学関数
# 三角関数
x = np.array([0, np.pi/6, np.pi/4, np.pi/3, np.pi/2])
print(np.sin(x)) # [0. 0.5 0.70710678 0.8660254 1. ]
print(np.cos(x)) # [1. 0.8660254 0.70710678 0.5 0. ]
# 指数・対数関数
y = np.array([1, 2, 3])
print(np.exp(y)) # [ 2.71828183 7.3890561 20.08553692]
print(np.log(y)) # [0. 0.69314718 1.09861229]
# 累乗・平方根
print(np.power(2, y)) # [2 4 8]
print(np.sqrt(y)) # [1. 1.41421356 1.73205081]
統計関数
data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
# 基本統計量
print(np.mean(data)) # 5.5
print(np.median(data)) # 5.5
print(np.std(data)) # 2.8722813232690143
print(np.var(data)) # 8.25
# 最大値・最小値とその位置
print(np.max(data)) # 10
print(np.min(data)) # 1
print(np.argmax(data)) # 9
print(np.argmin(data)) # 0
# 分位数
print(np.percentile(data, 25)) # 3.25
print(np.percentile(data, 75)) # 7.75
random モジュールによる乱数生成
import random
基本的な乱数生成
# 0以上1未満の乱数
print(random.random()) # 例: 0.6394267984578837
# 指定範囲の整数
print(random.randint(1, 10)) # 1から10までの整数
# 指定範囲の浮動小数点数
print(random.uniform(1.0, 10.0)) # 1.0から10.0までの実数
# リストからランダム選択
colors = ['red', 'blue', 'green', 'yellow']
print(random.choice(colors)) # 例: 'blue'
# 複数選択（重複あり）
print(random.choices(colors, k=3)) # 例: ['red', 'red', 'blue']
# 複数選択（重複なし）
print(random.sample(colors, 2)) # 例: ['green', 'red']
シャッフルと確率分布
# リストのシャッフル
numbers = [1, 2, 3, 4, 5]
random.shuffle(numbers)
print(numbers) # 例: [3, 1, 5, 2, 4]
# 正規分布
print(random.gauss(0, 1)) # 平均0、標準偏差1の正規分布
print(random.normalvariate(100, 15)) # 平均100、標準偏差15
# その他の分布
print(random.expovariate(1.5)) # 指数分布
print(random.betavariate(2, 5)) # ベータ分布
シード設定による再現性
# シードを設定して再現可能な乱数
random.seed(42)
print(random.random()) # 常に同じ値: 0.6394267984578837
random.seed(42)
print(random.random()) # 再び同じ値: 0.6394267984578837
fractions で分数計算
from fractions import Fraction
分数の基本操作
# 分数の作成
f1 = Fraction(3, 4) # 3/4
f2 = Fraction(1, 2) # 1/2
f3 = Fraction('0.25') # 文字列から: 1/4
f4 = Fraction(0.5) # 浮動小数点から: 1/2
print(f1) # 3/4
print(f3) # 1/4
# 分数の演算
print(f1 + f2) # 5/4
print(f1 - f2) # 1/4
print(f1 * f2) # 3/8
print(f1 / f2) # 3/2
# 分数と整数の演算
print(f1 + 1) # 7/4
print(f1 * 3) # 9/4
実用的な分数計算
# 料理のレシピ計算
recipe_flour = Fraction(2, 3) # 小麦粉 2/3カップ
recipe_sugar = Fraction(1, 4) # 砂糖 1/4カップ
# 1.5倍にする
scale = Fraction(3, 2)
print(f"小麦粉: {recipe_flour * scale}カップ") # 1カップ
print(f"砂糖: {recipe_sugar * scale}カップ") # 3/8カップ
# 約分と変換
f = Fraction(6, 8)
print(f) # 3/4（自動で約分）
print(float(f)) # 0.75
print(f.numerator) # 分子: 3
print(f.denominator) # 分母: 4
decimal で高精度小数計算
from decimal import Decimal, getcontext
高精度計算
# 浮動小数点の精度問題
print(0.1 + 0.2) # 0.30000000000000004
# Decimalで正確な計算
a = Decimal('0.1')
b = Decimal('0.2')
print(a + b) # 0.3
# 精度の設定
getcontext().prec = 50 # 50桁の精度
c = Decimal(1) / Decimal(3)
print(c) # 0.33333333333333333333333333333333333333333333333333
金融計算での活用
# 金融計算例
price = Decimal('19.99')
tax_rate = Decimal('0.08')
quantity = 3
subtotal = price * quantity
tax = subtotal * tax_rate
total = subtotal + tax
print(f"小計: ${subtotal}") # 小計: $59.97
print(f"税額: ${tax}") # 税額: $4.7976
print(f"合計: ${total}") # 合計: $64.7676
# 四捨五入
total_rounded = total.quantize(Decimal('0.01'))
print(f"合計（四捨五入）: ${total_rounded}") # 合計（四捨五入）: $64.77
組み合わせ・順列（itertools）
import itertools
組み合わせと順列
items = ['A', 'B', 'C', 'D']
# 組み合わせ（順序を考慮しない）
print("2つの組み合わせ:")
for combo in itertools.combinations(items, 2):
print(combo)
# ('A', 'B')
# ('A', 'C')
# ('A', 'D')
# ('B', 'C')
# ('B', 'D')
# ('C', 'D')
# 重複組み合わせ
print("\n重複を許す2つの組み合わせ:")
for combo in itertools.combinations_with_replacement(items[:2], 2):
print(combo)
# ('A', 'A')
# ('A', 'B')
# ('B', 'B')
# 順列（順序を考慮）
print("\n2つの順列:")
for perm in itertools.permutations(items, 2):
print(perm)
# ('A', 'B')
# ('A', 'C')
# ('A', 'D')
# ('B', 'A')
# ('B', 'C')
# ('B', 'D')
# ('C', 'A')
# ('C', 'B')
# ('C', 'D')
# ('D', 'A')
# ('D', 'B')
# ('D', 'C')
直積と繰り返し
# 直積
colors = ['red', 'blue']
sizes = ['S', 'M', 'L']
print("商品の組み合わせ:")
for product in itertools.product(colors, sizes):
print(product)
# ('red', 'S')
# ('red', 'M')
# ('red', 'L')
# ('blue', 'S')
# ('blue', 'M')
# ('blue', 'L')
# 数字の組み合わせ計算
def count_combinations(n, r):
"""組み合わせの数を計算"""
return len(list(itertools.combinations(range(n), r)))
print(f"5つから3つ選ぶ組み合わせ: {count_combinations(5, 3)}個") # 10個
複素数の計算
# 複素数の作成
z1 = 3 + 4j
z2 = 1 - 2j
# 基本演算
print(z1 + z2) # (4+2j)
print(z1 * z2) # (11-2j)
# 複素数の性質
print(f"実部: {z1.real}") # 3.0
print(f"虚部: {z1.imag}") # 4.0
print(f"絶対値: {abs(z1)}") # 5.0
print(f"共役: {z1.conjugate()}") # (3-4j)
# cmathモジュールで複素数関数
import cmath
print(f"位相: {cmath.phase(z1)}") # 0.9272952180016122
print(f"極形式: {cmath.polar(z1)}") # (5.0, 0.9272952180016122)
print(f"指数形式: {cmath.exp(1j * cmath.pi)}") # (-1+1.2246467991473532e-16j)
SciPyによる高度な数学処理
import scipy
from scipy import optimize, integrate, interpolate, special
最適化問題
# 関数の最小値を求める
def objective_function(x):
return x**2 + 2*x + 1
result = optimize.minimize_scalar(objective_function)
print(f"最小値: {result.fun} (x = {result.x})")
# 方程式の解を求める
def equation(x):
return x**3 - 2*x - 5
solution = optimize.fsolve(equation, 2)
print(f"方程式の解: x = {solution[0]}")
数値積分
# 定積分の計算
def integrand(x):
return x**2
result, error = integrate.quad(integrand, 0, 1)
print(f"∫₀¹ x² dx = {result} (誤差: {error})")
# 重積分
def integrand_2d(y, x):
return x * y
result_2d, error_2d = integrate.dblquad(integrand_2d, 0, 1, 0, 1)
print(f"二重積分の結果: {result_2d}")
特殊関数
# ガンマ関数
print(f"Γ(5) = {special.gamma(5)}") # 4! = 24
# ベッセル関数
x = np.linspace(0, 10, 100)
j0 = special.jv(0, x) # 第1種ベッセル関数
# 誤差関数
print(f"erf(1) = {special.erf(1)}")
SymPyによる記号計算
import sympy as sp
記号の定義と基本操作
# 記号の定義
x, y = sp.symbols('x y')
# 式の定義
expr = x**2 + 2*x + 1
print(expr) # x**2 + 2*x + 1
# 式の展開
expanded = sp.expand((x + 1)**3)
print(expanded) # x**3 + 3*x**2 + 3*x + 1
# 因数分解
factored = sp.factor(x**2 - 4)
print(factored) # (x - 2)*(x + 2)
微分と積分
# 微分
f = x**3 + 2*x**2 + x
df_dx = sp.diff(f, x)
print(f"f'(x) = {df_dx}") # 3*x**2 + 4*x + 1
# 積分
integral = sp.integrate(f, x)
print(f"∫f(x)dx = {integral}") # x**4/4 + 2*x**3/3 + x**2/2
# 定積分
definite_integral = sp.integrate(x**2, (x, 0, 1))
print(f"∫₀¹ x² dx = {definite_integral}") # 1/3
方程式の解
# 1次方程式
eq1 = sp.Eq(2*x + 3, 7)
sol1 = sp.solve(eq1, x)
print(f"2x + 3 = 7 の解: x = {sol1}") # [2]
# 2次方程式
eq2 = x**2 - 5*x + 6
sol2 = sp.solve(eq2, x)
print(f"x² - 5x + 6 = 0 の解: x = {sol2}") # [2, 3]
# 連立方程式
eq3 = [x + y - 5, 2*x - y - 1]
sol3 = sp.solve(eq3, [x, y])
print(f"連立方程式の解: {sol3}") # {x: 2, y: 3}
実用的な数学処理の例
フィボナッチ数列
def fibonacci_sequence(n):
"""フィボナッチ数列の最初のn項を返す"""
fib = [0, 1]
for i in range(2, n):
fib.append(fib[i-1] + fib[i-2])
return fib[:n]
# 黄金比の近似
fib_nums = fibonacci_sequence(20)
ratios = [fib_nums[i]/fib_nums[i-1] for i in range(2, len(fib_nums))]
print("フィボナッチ数列の比（黄金比に収束）:")
for i, ratio in enumerate(ratios[-5:], start=len(ratios)-4):
print(f"F({i+2})/F({i+1}) = {ratio:.6f}")
golden_ratio = (1 + math.sqrt(5)) / 2
print(f"黄金比の理論値: {golden_ratio:.6f}")
素数判定と生成
def is_prime(n):
"""素数判定"""
if n < 2:
return False
for i in range(2, int(math.sqrt(n)) + 1):
if n % i == 0:
return False
return True
def sieve_of_eratosthenes(limit):
"""エラトステネスの篩による素数生成"""
is_prime_arr = [True] * (limit + 1)
is_prime_arr[0] = is_prime_arr[1] = False
for i in range(2, int(math.sqrt(limit)) + 1):
if is_prime_arr[i]:
for j in range(i*i, limit + 1, i):
is_prime_arr[j] = False
return [i for i in range(2, limit + 1) if is_prime_arr[i]]
# 100以下の素数
primes = sieve_of_eratosthenes(100)
print(f"100以下の素数: {primes}")
print(f"100以下の素数の個数: {len(primes)}")
モンテカルロ法で円周率を求める
import random
def estimate_pi(n_samples):
"""モンテカルロ法で円周率を推定"""
inside_circle = 0
for _ in range(n_samples):
x = random.uniform(-1, 1)
y = random.uniform(-1, 1)
if x*x + y*y <= 1:
inside_circle += 1
return 4 * inside_circle / n_samples
# 円周率の推定
for n in [1000, 10000, 100000, 1000000]:
pi_estimate = estimate_pi(n)
error = abs(pi_estimate - math.pi)
print(f"サンプル数 {n:7d}: π ≈ {pi_estimate:.6f} (誤差: {error:.6f})")
数値微分
def numerical_derivative(f, x, h=1e-8):
"""数値微分（中央差分）"""
return (f(x + h) - f(x - h)) / (2 * h)
# 関数とその導関数
def f(x):
return x**3 + 2*x**2 + x
def f_prime_exact(x):
return 3*x**2 + 4*x + 1
# 数値微分の精度確認
x = 2.0
numerical = numerical_derivative(f, x)
exact = f_prime_exact(x)
print(f"x = {x} での微分:")
print(f"数値微分: {numerical:.8f}")
print(f"解析解: {exact:.8f}")
print(f"誤差: {abs(numerical - exact):.2e}")
ニュートン法による方程式の解
def newton_method(f, df, x0, tolerance=1e-10, max_iterations=100):
"""ニュートン法による方程式の解"""
x = x0
for i in range(max_iterations):
fx = f(x)
if abs(fx) < tolerance:
return x, i
dfx = df(x)
if abs(dfx) < tolerance:
raise ValueError("導関数がゼロに近すぎます")
x = x - fx / dfx
raise ValueError("収束しませんでした")
# √2 を求める（x² - 2 = 0 の解）
def f(x):
return x**2 - 2
def df(x):
return 2*x
root, iterations = newton_method(f, df, 1.0)
print(f"√2 ≈ {root:.10f} ({iterations}回の反復)")
print(f"理論値: {math.sqrt(2):.10f}")
print(f"誤差: {abs(root - math.sqrt(2)):.2e}")
行列と線形代数の応用
連立一次方程式の解
# 連立一次方程式 Ax = b を解く
# 2x + 3y = 7
# x - y = 1
A = np.array([[2, 3], [1, -1]])
b = np.array([7, 1])
# 解を求める
x = np.linalg.solve(A, b)
print(f"解: x = {x[0]}, y = {x[1]}")
# 検算
print(f"検算: Ax = {A @ x}")
print(f"b = {b}")
主成分分析（PCA）の基本
# サンプルデータ
np.random.seed(42)
data = np.random.randn(100, 2)
data[:, 1] = data[:, 0] * 2 + np.random.randn(100) * 0.5
# 共分散行列
cov_matrix = np.cov(data.T)
print("共分散行列:")
print(cov_matrix)
# 固有値と固有ベクトル
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
print(f"固有値: {eigenvalues}")
print(f"固有ベクトル:")
print(eigenvectors)
# 寄与率
contribution_rate = eigenvalues / eigenvalues.sum()
print(f"寄与率: {contribution_rate}")
最適化問題の実例
線形計画法
from scipy.optimize import linprog
# 最大化問題: 3x + 2y を最大化
# 制約条件: x + y ≤ 4, 2x + y ≤ 5, x ≥ 0, y ≥ 0
# linprogは最小化なので、符号を反転
c = [-3, -2] # 目的関数の係数（符号反転）
A_ub = [[1, 1], [2, 1]] # 不等式制約の係数行列
b_ub = [4, 5] # 不等式制約の右辺
bounds = [(0, None), (0, None)] # 変数の範囲
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')
if result.success:
print(f"最適解: x = {result.x[0]:.2f}, y = {result.x[1]:.2f}")
print(f"最大値: {-result.fun:.2f}")
else:
print("最適化に失敗しました")
非線形最適化
from scipy.optimize import minimize
# 目的関数: (x-3)² + (y-2)² を最小化
def objective(vars):
x, y = vars
return (x - 3)**2 + (y - 2)**2
# 制約条件: x + y ≥ 1
def constraint(vars):
x, y = vars
return x + y - 1
# 制約を定義
con = {'type': 'ineq', 'fun': constraint}
# 初期値
x0 = [0, 0]
# 最適化実行
result = minimize(objective, x0, method='SLSQP', constraints=con)
if result.success:
print(f"最適解: x = {result.x[0]:.4f}, y = {result.x[1]:.4f}")
print(f"最小値: {result.fun:.4f}")
数値解析の高度なテクニック
ルンゲ・クッタ法による微分方程式の数値解
def runge_kutta_4(f, y0, t0, t1, n):
"""4次ルンゲ・クッタ法"""
h = (t1 - t0) / n
t = t0
y = y0
t_values = [t0]
y_values = [y0]
for i in range(n):
k1 = h * f(t, y)
k2 = h * f(t + h/2, y + k1/2)
k3 = h * f(t + h/2, y + k2/2)
k4 = h * f(t + h, y + k3)
y = y + (k1 + 2*k2 + 2*k3 + k4) / 6
t = t + h
t_values.append(t)
y_values.append(y)
return np.array(t_values), np.array(y_values)
# 微分方程式 dy/dt = -2y + 1, y(0) = 0 の解
def dydt(t, y):
return -2*y + 1
t_vals, y_vals = runge_kutta_4(dydt, 0, 0, 2, 100)
# 解析解との比較
def analytical_solution(t):
return 0.5 * (1 - np.exp(-2*t))
y_analytical = analytical_solution(t_vals)
error = np.abs(y_vals - y_analytical)
print(f"t=2での数値解: {y_vals[-1]:.6f}")
print(f"t=2での解析解: {y_analytical[-1]:.6f}")
print(f"最大誤差: {np.max(error):.2e}")
フーリエ変換
# 信号の生成
t = np.linspace(0, 1, 1000)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 10 * t)
# ノイズを追加
noise = 0.2 * np.random.randn(len(t))
noisy_signal = signal + noise
# フーリエ変換
fft_result = np.fft.fft(noisy_signal)
frequencies = np.fft.fftfreq(len(t), t[1] - t[0])
# パワースペクトルの計算
power_spectrum = np.abs(fft_result)**2
# ピーク周波数の検出
positive_freq_idx = frequencies > 0
peak_indices = np.argsort(power_spectrum[positive_freq_idx])[-2:]
peak_frequencies = frequencies[positive_freq_idx][peak_indices]
print(f"検出された周波数: {sorted(peak_frequencies):.1f} Hz")
print("理論値: [5.0, 10.0] Hz")
統計学と確率論の応用
中心極限定理の確認
def demonstrate_central_limit_theorem(population_dist, sample_size, num_samples):
"""中心極限定理のデモンストレーション"""
sample_means = []
for _ in range(num_samples):
sample = np.random.choice(population_dist, sample_size)
sample_means.append(np.mean(sample))
return np.array(sample_means)
# 一様分布からのサンプリング
population = np.random.uniform(0, 10, 10000)
# 異なるサンプルサイズでの検証
sample_sizes = [1, 5, 30, 100]
num_samples = 1000
print("中心極限定理の確認:")
print(f"母集団の平均: {np.mean(population):.2f}")
print(f"母集団の標準偏差: {np.std(population):.2f}")
print()
for n in sample_sizes:
sample_means = demonstrate_central_limit_theorem(population, n, num_samples)
theoretical_std = np.std(population) / np.sqrt(n)
print(f"サンプルサイズ {n}:")
print(f" 標本平均の平均: {np.mean(sample_means):.3f}")
print(f" 標本平均の標準偏差: {np.std(sample_means):.3f}")
print(f" 理論的標準偏差: {theoretical_std:.3f}")
print()
ベイズ統計の基本
def bayes_theorem(prior, likelihood, evidence):
"""ベイズの定理"""
return (likelihood * prior) / evidence
# 例: 病気の診断
# 病気の事前確率: 1%
prior_disease = 0.01
prior_healthy = 0.99
# 検査の精度: 病気の人が陽性になる確率95%, 健康な人が陽性になる確率5%
likelihood_positive_given_disease = 0.95
likelihood_positive_given_healthy = 0.05
# 陽性になる全体的な確率
evidence_positive = (likelihood_positive_given_disease * prior_disease +
likelihood_positive_given_healthy * prior_healthy)
# 陽性の場合に実際に病気である確率
posterior_disease = bayes_theorem(prior_disease,
likelihood_positive_given_disease,
evidence_positive)
print("医療診断のベイズ統計例:")
print(f"病気の事前確率: {prior_disease:.1%}")
print(f"検査陽性の場合の病気確率: {posterior_disease:.1%}")
機械学習の数学的基礎
線形回帰の実装
def linear_regression_normal_equation(X, y):
"""正規方程式による線形回帰"""
# X に切片項を追加
X_with_intercept = np.column_stack([np.ones(len(X)), X])
# 正規方程式: θ = (X^T X)^(-1) X^T y
theta = np.linalg.inv(X_with_intercept.T @ X_with_intercept) @ X_with_intercept.T @ y
return theta
# サンプルデータ
np.random.seed(42)
X = np.linspace(0, 10, 100)
y = 2 * X + 1 + np.random.randn(100) * 0.5
# 線形回帰の実行
theta = linear_regression_normal_equation(X, y)
intercept, slope = theta
print(f"推定された切片: {intercept:.3f}")
print(f"推定された傾き: {slope:.3f}")
print("理論値: 切片=1.0, 傾き=2.0")
# 予測
y_pred = intercept + slope * X
# 決定係数（R²）の計算
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
r_squared = 1 - (ss_res / ss_tot)
print(f"決定係数 R²: {r_squared:.4f}")
勾配降下法の実装
def gradient_descent(X, y, learning_rate=0.01, num_iterations=1000):
"""勾配降下法による線形回帰"""
m = len(X)
X_with_intercept = np.column_stack([np.ones(m), X])
theta = np.zeros(2) # 初期値
cost_history = []
for i in range(num_iterations):
# 予測
h = X_with_intercept @ theta
# コスト計算
cost = (1 / (2 * m)) * np.sum((h - y) ** 2)
cost_history.append(cost)
# 勾配計算
gradient = (1 / m) * X_with_intercept.T @ (h - y)
# パラメータ更新
theta = theta - learning_rate * gradient
return theta, cost_history
# 勾配降下法の実行
theta_gd, costs = gradient_descent(X, y, learning_rate=0.001, num_iterations=1000)
print(f"勾配降下法による推定:")
print(f"切片: {theta_gd[0]:.3f}, 傾き: {theta_gd[1]:.3f}")
print(f"最終コスト: {costs[-1]:.4f}")
まとめ：Pythonで数学の幅が広がる！
Pythonには、数値処理に必要な機能がほぼすべて揃っています。この記事で学んだ内容をまとめると：
標準ライブラリによる基本処理
- math：基本的な数学関数と定数
- statistics：統計計算
- random：乱数生成
- fractions：分数計算
- decimal：高精度小数計算
- itertools：組み合わせ・順列
外部ライブラリによる高度な処理
- NumPy：配列演算と線形代数
- SciPy：科学計算全般
- SymPy：記号計算
- matplotlib：グラフ描画（可視化）
実用的な応用分野
- 統計解析とデータサイエンス
- 機械学習の数学的基礎
- 数値解析と最適化
- 信号処理とフーリエ解析
- 確率論とベイズ統計
パフォーマンスのポイント
- NumPyを使ったベクトル化計算
- 適切なデータ型の選択
- メモリ効率を考慮したアルゴリズム
- 必要に応じた並列処理
コメント