[Python/Scipy] curve_fitの使い方 – カーブフィッティング手法
curve_fit
は、PythonのSciPyライブラリに含まれる関数で、非線形のカーブフィッティングを行うために使用されます。
与えられたデータに対して、指定した関数のパラメータを最適化し、最も適合する曲線を見つけます。
基本的な使い方は、まずフィッティングしたい関数を定義し、次にcurve_fit
にその関数、データのx値、y値を渡します。
返り値として、最適なパラメータとその共分散行列が得られます。
curve_fitとは?基本的な概要
curve_fit
は、PythonのSciPyライブラリに含まれる関数で、データに対して最適な曲線をフィッティングするためのツールです。
この関数は、与えられたデータポイントに基づいて、指定したモデル関数のパラメータを最適化します。
フィッティングの目的は、データのトレンドを把握したり、将来の値を予測したりすることです。
curve_fit
は、最小二乗法を用いてパラメータを推定し、データの誤差を最小化します。
これにより、実際のデータに最も適した曲線を見つけることができます。
例えば、線形関数や非線形関数、指数関数、ガウス関数など、さまざまなモデルを使用してフィッティングを行うことが可能です。
この機能は、科学技術計算やデータ分析の分野で広く利用されており、実験データの解析や機械学習の前処理など、多岐にわたる応用が期待されています。
curve_fitの基本的な使い方
フィッティング関数の定義
フィッティングを行うためには、まずモデルとなる関数を定義する必要があります。
例えば、線形関数や非線形関数を定義します。
以下は、線形関数の例です。
def linear_function(x, a, b):
return a * x + b
この関数は、傾き\(a\)と切片\(b\)を持つ直線を表します。
データの準備
次に、フィッティングに使用するデータを準備します。
データは通常、x座標とy座標のペアとして用意されます。
以下は、サンプルデータの生成例です。
import numpy as np
# サンプルデータの生成
x_data = np.array([1, 2, 3, 4, 5])
y_data = np.array([2.2, 2.8, 3.6, 4.5, 5.1])
curve_fitの呼び出し
curve_fit関数
を使用して、データにフィッティングを行います。
以下のコードでは、先ほど定義した線形関数を用いてフィッティングを行います。
from scipy.optimize import curve_fit
# curve_fitの呼び出し
popt, pcov = curve_fit(linear_function, x_data, y_data)
ここで、popt
には最適化されたパラメータが、pcov
には共分散行列が格納されます。
最適化されたパラメータの取得
フィッティングの結果、最適化されたパラメータを取得することができます。
以下のコードで、傾きと切片を表示します。
a_opt, b_opt = popt
print(f"最適化された傾き: {a_opt}, 最適化された切片: {b_opt}")
出力結果は以下のようになります。
最適化された傾き: 0.7499999995497489, 最適化された切片: 1.3900000016509204
共分散行列の取得と解釈
共分散行列は、フィッティングパラメータの不確かさを示します。
pcov
を使って、各パラメータの標準誤差を計算することができます。
以下のコードで、標準誤差を表示します。
import np
# 標準誤差の計算
perr = np.sqrt(np.diag(pcov))
print(f"傾きの標準誤差: {perr[0]}, 切片の標準誤差: {perr[1]}")
出力結果は以下のようになります。
傾きの標準誤差: 0.030000000578937647, 切片の標準誤差: 0.09949874366827378
このように、共分散行列を用いることで、フィッティングパラメータの信頼性を評価することができます。
実際の例:線形フィッティング
線形関数の定義
線形フィッティングを行うために、まず線形関数を定義します。
この関数は、傾きと切片を持つ直線を表します。
以下のように定義します。
def linear_function(x, a, b):
return a * x + b
サンプルデータの生成
次に、フィッティングに使用するサンプルデータを生成します。
ここでは、ノイズを加えた線形データを作成します。
import numpy as np
# サンプルデータの生成
np.random.seed(0) # 再現性のための乱数シード
x_data = np.linspace(0, 10, 20) # 0から10までの20点
true_a = 2.0 # 真の傾き
true_b = 1.0 # 真の切片
y_data = true_a * x_data + true_b + np.random.normal(0, 1, size=x_data.shape) # ノイズを加えたデータ
curve_fitを使ったフィッティング
生成したデータに対して、curve_fit
を使用してフィッティングを行います。
以下のコードで、最適化されたパラメータを取得します。
from scipy.optimize import curve_fit
# curve_fitの呼び出し
popt, pcov = curve_fit(linear_function, x_data, y_data)
# 最適化されたパラメータの取得
a_opt, b_opt = popt
print(f"最適化された傾き: {a_opt}, 最適化された切片: {b_opt}")
出力結果は以下のようになります。
最適化された傾き: 1.8862774176485506, 最適化された切片: 2.137947497608783
フィッティング結果の可視化
最後に、フィッティング結果を可視化します。
Matplotlibを使用して、元のデータとフィッティングした直線をプロットします。
import matplotlib.pyplot as plt
# フィッティング結果の可視化
plt.scatter(x_data, y_data, label='データ', color='blue')
plt.plot(x_data, linear_function(x_data, a_opt, b_opt), label='フィッティング直線', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('線形フィッティングの例')
plt.legend()
plt.grid()
plt.show()
このコードを実行すると、元のデータポイントとフィッティングした直線が表示されます。
これにより、フィッティングの結果を視覚的に確認することができます。
実際の例:非線形フィッティング
非線形関数の定義
非線形フィッティングを行うために、まず非線形関数を定義します。
ここでは、指数関数を用いたモデルを考えます。
以下のように定義します。
def exponential_function(x, a, b):
return a * np.exp(b * x)
この関数は、パラメータ\(a\)と\(b\)を持つ指数関数を表します。
サンプルデータの生成
次に、フィッティングに使用するサンプルデータを生成します。
ここでは、真のパラメータを持つ指数関数にノイズを加えたデータを作成します。
import numpy as np
# サンプルデータの生成
np.random.seed(0) # 再現性のための乱数シード
x_data = np.linspace(0, 5, 50) # 0から5までの50点
true_a = 2.0 # 真のパラメータa
true_b = 1.5 # 真のパラメータb
y_data = true_a * np.exp(true_b * x_data) + np.random.normal(0, 0.5, size=x_data.shape) # ノイズを加えたデータ
curve_fitを使ったフィッティング
生成したデータに対して、curve_fit
を使用してフィッティングを行います。
以下のコードで、最適化されたパラメータを取得します。
from scipy.optimize import curve_fit
# curve_fitの呼び出し
popt, pcov = curve_fit(exponential_function, x_data, y_data)
# 最適化されたパラメータの取得
a_opt, b_opt = popt
print(f"最適化されたパラメータa: {a_opt}, 最適化されたパラメータb: {b_opt}")
出力結果は以下のようになります。
最適化されたパラメータa: 1.9998422329660035, 最適化されたパラメータb: 1.4999963109691483
フィッティング結果の可視化
最後に、フィッティング結果を可視化します。
Matplotlibを使用して、元のデータとフィッティングした指数関数をプロットします。
import matplotlib.pyplot as plt
# フィッティング結果の可視化
plt.scatter(x_data, y_data, label='データ', color='blue')
plt.plot(x_data, exponential_function(x_data, a_opt, b_opt), label='フィッティング曲線', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('非線形フィッティングの例')
plt.legend()
plt.grid()
plt.show()
このコードを実行すると、元のデータポイントとフィッティングした指数関数が表示されます。
これにより、非線形フィッティングの結果を視覚的に確認することができます。
curve_fitのオプションとパラメータ
curve_fit関数
には、フィッティングの精度や効率を向上させるためのさまざまなオプションやパラメータがあります。
以下に主要なオプションを説明します。
初期パラメータの設定
フィッティングを行う際、初期パラメータを設定することができます。
これにより、最適化プロセスがより効率的に行われる場合があります。
初期パラメータは、p0
引数を使用して指定します。
以下の例では、初期パラメータを設定しています。
# 初期パラメータの設定
initial_params = [1.0, 1.0] # a=1.0, b=1.0
popt, pcov = curve_fit(exponential_function, x_data, y_data, p0=initial_params)
境界条件の設定
フィッティングパラメータに対して境界条件を設定することも可能です。
これにより、パラメータが特定の範囲内に収まるように制約をかけることができます。
境界条件は、bounds
引数を使用して指定します。
以下の例では、パラメータの範囲を設定しています。
# 境界条件の設定
bounds = (0, [3.0, 2.0]) # aは0以上、bは2.0以下
popt, pcov = curve_fit(exponential_function, x_data, y_data, bounds=bounds)
絶対誤差と相対誤差の設定
データの誤差を考慮するために、絶対誤差や相対誤差を設定することができます。
これにより、フィッティングの精度を向上させることができます。
誤差は、sigma
引数を使用して指定します。
以下の例では、絶対誤差を設定しています。
# 絶対誤差の設定
absolute_error = 0.5 # 各データポイントの誤差
popt, pcov = curve_fit(exponential_function, x_data, y_data, sigma=absolute_error)
最大反復回数の設定
最適化プロセスの最大反復回数を設定することも可能です。
これにより、フィッティングが収束しない場合に、無限ループを防ぐことができます。
最大反復回数は、maxfev
引数を使用して指定します。
以下の例では、最大反復回数を設定しています。
# 最大反復回数の設定
max_iterations = 1000
popt, pcov = curve_fit(exponential_function, x_data, y_data, maxfev=max_iterations)
これらのオプションを適切に設定することで、curve_fit
をより効果的に活用し、フィッティングの精度を向上させることができます。
フィッティングの精度を向上させる方法
フィッティングの精度を向上させるためには、いくつかの方法があります。
以下に、具体的なアプローチを紹介します。
初期パラメータの選び方
初期パラメータは、最適化プロセスの収束に大きな影響を与えます。
適切な初期パラメータを選ぶためには、以下のポイントを考慮します。
- データの可視化: データをプロットして、トレンドを視覚的に確認します。
これにより、初期パラメータの推定が容易になります。
- 経験則の利用: 過去のデータや類似の問題から得た経験則を基に、初期パラメータを設定します。
- 試行錯誤: 異なる初期パラメータを試して、最も良い結果を得られるものを選びます。
データの前処理
データの前処理は、フィッティングの精度を向上させるために重要です。
以下の手法を考慮します。
- 外れ値の除去: 外れ値がフィッティング結果に悪影響を与えることがあるため、外れ値を特定して除去します。
- スムージング: ノイズを軽減するために、データをスムージングする手法(移動平均やローパスフィルタなど)を適用します。
- 正規化: データのスケールを揃えるために、正規化を行います。
これにより、フィッティングの安定性が向上します。
ノイズの影響を軽減する方法
ノイズの影響を軽減するためには、以下の方法を検討します。
- データ収集の改善: 測定機器の精度を向上させたり、測定条件を最適化することで、ノイズを減少させます。
- 重み付けフィッティング:
curve_fit
のsigma
引数を使用して、データポイントに重みを付けることで、ノイズの影響を考慮します。 - フィルタリング: データに対してフィルタリングを行い、ノイズを除去します。
例えば、バターワースフィルタやカラーマンフィルタなどが利用できます。
フィッティング結果の評価方法
フィッティング結果を評価することで、精度を確認し、必要に応じて改善策を講じることができます。
以下の評価方法を考慮します。
- 残差分析: フィッティング結果と実データの差(残差)を分析し、パターンや外れ値を特定します。
残差がランダムであれば、フィッティングが適切であると考えられます。
- 決定係数(R²): フィッティングの良さを示す指標として、決定係数を計算します。
値が1に近いほど、フィッティングが良好であることを示します。
- 交差検証: データをトレーニングセットとテストセットに分け、フィッティングの精度を評価します。
これにより、過学習を防ぐことができます。
これらの方法を組み合わせて実施することで、フィッティングの精度を向上させることが可能です。
応用例:複数の変数を持つ関数のフィッティング
複数の変数を持つ関数のフィッティングは、実際のデータ分析において非常に重要です。
ここでは、2つの変数を持つ関数のフィッティングの例を示します。
複数変数の関数定義
まず、2つの変数を持つ関数を定義します。
ここでは、2次元の多項式関数を例にします。
以下のように定義します。
def polynomial_function(X, a, b, c):
x, y = X
return a * x**2 + b * y**2 + c
この関数は、\(x\)と\(y\)の2つの変数に対して、係数\(a\)、\(b\)、および定数項\(c\)を持つ2次の多項式を表します。
データの準備
次に、フィッティングに使用するデータを準備します。
ここでは、ノイズを加えたサンプルデータを生成します。
import numpy as np
# サンプルデータの生成
np.random.seed(0) # 再現性のための乱数シード
x_data = np.linspace(-5, 5, 20) # xの範囲
y_data = np.linspace(-5, 5, 20) # yの範囲
x_data, y_data = np.meshgrid(x_data, y_data) # グリッドデータの生成
# 真のパラメータ
true_a = 1.0
true_b = 2.0
true_c = 3.0
# ノイズを加えたデータの生成
z_data = true_a * x_data**2 + true_b * y_data**2 + true_c + np.random.normal(0, 5, size=x_data.shape)
curve_fitを使ったフィッティング
生成したデータに対して、curve_fit
を使用してフィッティングを行います。
curve_fit
は、2次元のデータに対しても使用できますが、引数として1次元の配列を必要とするため、データをフラットにする必要があります。
from scipy.optimize import curve_fit
# データをフラットにする
x_flat = x_data.flatten()
y_flat = y_data.flatten()
z_flat = z_data.flatten()
# curve_fitの呼び出し
popt, pcov = curve_fit(polynomial_function, (x_flat, y_flat), z_flat)
# 最適化されたパラメータの取得
a_opt, b_opt, c_opt = popt
print(f"最適化されたパラメータ a: {a_opt}, b: {b_opt}, c: {c_opt}")
出力結果は以下のようになります。
最適化されたパラメータ a: 0.9868686176365955, b: 2.007189926308829, c: 2.9119639603638277
フィッティング結果の解釈
フィッティング結果を解釈するためには、最適化されたパラメータを確認し、元のデータとフィッティングした関数を比較します。
最適化されたパラメータが真のパラメータに近い場合、フィッティングが成功したと考えられます。
また、フィッティング結果を可視化することで、フィッティングの良さを確認できます。
以下のコードで、元のデータとフィッティングした関数をプロットします。
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# フィッティング結果の可視化
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x_flat, y_flat, z_flat, label='データ', color='blue')
# フィッティングした関数の計算
z_fit = polynomial_function((x_flat, y_flat), a_opt, b_opt, c_opt).reshape(x_data.shape)
# フィッティング曲面のプロット
ax.plot_surface(x_data, y_data, z_fit, color='red', alpha=0.5, label='フィッティング曲面')
ax.set_xlabel('X軸')
ax.set_ylabel('Y軸')
ax.set_zlabel('Z軸')
ax.set_title('複数変数のフィッティング結果')
plt.legend()
plt.show()
このプロットにより、元のデータとフィッティングした関数の関係を視覚的に確認することができます。
フィッティングの精度を評価し、必要に応じてモデルを改善することが可能です。
応用例:指数関数フィッティング
指数関数フィッティングは、成長や減衰の過程をモデル化する際に非常に有用です。
ここでは、実データを用いて指数関数のフィッティングを行う例を示します。
指数関数の定義
まず、フィッティングに使用する指数関数を定義します。
以下のように、パラメータ\(a\)と\(b\)を持つ指数関数を定義します。
def exponential_function(x, a, b):
return a * np.exp(b * x)
この関数は、\(a\)が初期値、\(b\)が成長率を表します。
実データを使ったフィッティング
次に、実データを用いてフィッティングを行います。
ここでは、サンプルデータを生成し、それに対してフィッティングを行います。
実際のデータを使用する場合は、データを適切に読み込む必要があります。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# サンプルデータの生成
np.random.seed(0) # 再現性のための乱数シード
x_data = np.linspace(0, 5, 50) # xの範囲
true_a = 2.0 # 真のパラメータa
true_b = 1.5 # 真のパラメータb
y_data = true_a * np.exp(true_b * x_data) + np.random.normal(0, 0.5, size=x_data.shape) # ノイズを加えたデータ
# curve_fitを使ったフィッティング
popt, pcov = curve_fit(exponential_function, x_data, y_data)
# 最適化されたパラメータの取得
a_opt, b_opt = popt
print(f"最適化されたパラメータ a: {a_opt}, b: {b_opt}")
出力結果は以下のようになります。
最適化されたパラメータ a: 1.9998422329660035, b: 1.4999963109691483
フィッティング結果の可視化
最後に、フィッティング結果を可視化します。
元のデータとフィッティングした指数関数をプロットして、フィッティングの良さを確認します。
# フィッティング結果の可視化
plt.scatter(x_data, y_data, label='データ', color='blue')
plt.plot(x_data, exponential_function(x_data, a_opt, b_opt), label='フィッティング曲線', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('指数関数フィッティングの例')
plt.legend()
plt.grid()
plt.show()
このコードを実行すると、元のデータポイントとフィッティングした指数関数が表示されます。
これにより、フィッティングの結果を視覚的に確認することができ、モデルの適合度を評価することができます。
フィッティングの精度が高い場合、フィッティング曲線はデータポイントに非常に近い位置に表示されるはずです。
応用例:ガウス関数フィッティング
ガウス関数フィッティングは、データが正規分布に従う場合やピークを持つデータの解析に非常に有用です。
ここでは、ガウス関数を用いたフィッティングの例を示します。
ガウス関数の定義
まず、フィッティングに使用するガウス関数を定義します。
以下のように、平均値\( \mu \)と標準偏差\( \sigma \)を持つガウス関数を定義します。
def gaussian_function(x, amp, mu, sigma):
return amp * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2))
この関数は、振幅\( \text{amp} \)、平均値\( \mu \)、および標準偏差\( \sigma \)をパラメータとして持ちます。
実データを使ったフィッティング
次に、実データを用いてフィッティングを行います。
ここでは、サンプルデータを生成し、それに対してフィッティングを行います。
実際のデータを使用する場合は、データを適切に読み込む必要があります。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# サンプルデータの生成
np.random.seed(0) # 再現性のための乱数シード
x_data = np.linspace(-5, 5, 100) # xの範囲
true_amp = 1.0 # 真の振幅
true_mu = 0.0 # 真の平均値
true_sigma = 1.0 # 真の標準偏差
y_data = true_amp * np.exp(-((x_data - true_mu) ** 2) / (2 * true_sigma ** 2)) + np.random.normal(0, 0.05, size=x_data.shape) # ノイズを加えたデータ
# curve_fitを使ったフィッティング
popt, pcov = curve_fit(gaussian_function, x_data, y_data, p0=[1, 0, 1])
# 最適化されたパラメータの取得
amp_opt, mu_opt, sigma_opt = popt
print(f"最適化された振幅: {amp_opt}, 平均値: {mu_opt}, 標準偏差: {sigma_opt}")
出力結果は以下のようになります。
最適化された振幅: 0.9792553934510548, 平均値: -0.009608368753774257, 標準偏差: 0.9939791155630869
フィッティング結果の可視化
最後に、フィッティング結果を可視化します。
元のデータとフィッティングしたガウス関数をプロットして、フィッティングの良さを確認します。
# フィッティング結果の可視化
plt.scatter(x_data, y_data, label='データ', color='blue', s=10)
plt.plot(x_data, gaussian_function(x_data, amp_opt, mu_opt, sigma_opt), label='フィッティング曲線', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('ガウス関数フィッティングの例')
plt.legend()
plt.grid()
plt.show()
このコードを実行すると、元のデータポイントとフィッティングしたガウス関数が表示されます。
これにより、フィッティングの結果を視覚的に確認することができ、モデルの適合度を評価することができます。
フィッティングの精度が高い場合、フィッティング曲線はデータポイントに非常に近い位置に表示されるはずです。
まとめ
この記事では、PythonのSciPyライブラリにおけるcurve_fit
の使い方や、さまざまなフィッティング手法について詳しく解説しました。
特に、線形フィッティングや非線形フィッティング、複数変数を持つ関数のフィッティング、指数関数やガウス関数のフィッティングに関する具体的な例を通じて、実際のデータに対するフィッティングのプロセスを理解することができました。
これらの手法を活用して、実データの解析やモデルの構築に挑戦してみてください。